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INVESTIGATION OF OCCULTATION AND IMAGERY PROBLEMS 
FOR ORBITAL MISSIONS TO VENUS AND MARS 

By Richard N. Green 
Langley Research Center 

SUMMARY 

A computer program for the purpose of providing mission design data for planeto- 
centric elliptical orbits about either Mars or Venus has been developed. The times at 
which the satellite enters into and exits from the shadows of the Sun, Earth, and Canopus 
are computed. In addition, the positions in orbit that correspond to desirable lighting 
conditions on the surface for vertical photography are computed. A program option allows 
for these conditions to be examined at various times by permitting the planet to proceed 
in its orbit about the Sun. The mean heliocentric position of the planet, the angles in the 
Earth-Sun-planet triangle, and the associated distances are also included in the list of 
output parameters. The program is written in FORTRAN IV language for a digital com- 
puter and contains 27 subroutines of general use. A description of the program input and 
output, a complete FORTRAN listing of the program, and samples of the input and output 
are included. 


INTRODUCTION 

The computer program VAMOOS (Venus and Mars Orbital Occupation Simulator) 
originated for the purpose of studying orbital missions to Mars and Venus. Two impor- 
tant considerations to mission analysis are the occultations of celestial bodies and the 
available imagery. The length of time that the satellite spends in the shadows of the Sun, 
Earth, and Canopus is significant. When the satellite passes into the shadow of the Sun, 
the use of solar cells is prohibited. Many times this condition is undesirable or must be 
limited to short periods of darkness. Passage into the shadow of the Earth will result in 
the loss of communication between the satellite and Earth. Since the attitude reference 
system of the satellite is frequently alined with the Sun and Canopus, the occultation of 
either is of great interest. In addition to the importance of the length of time, the posi- 
tions in the orbit at which the satellite enters into and exits from the shadows are signifi- 
cant since they influence the operational sequence of events. The positions in the orbit 
of the satellite that correspond to desirable lighting conditions on the surface of the planet 
for vertical photography are also of interest. These positions, along with the orbital 



elements, lead to the definition of the surface which can be photographed. In addition, the 
resulting resolution is defined. If the mission objectives require that specific sections of 
the planet's surface be photographed, the orientation of the orbital plane of the satellite 
can be adjusted to satisfy the objectives. For these reasons VAMOOS was developed to 
supply the data necessary to define the occultations and the conditions at specific lighting 
angles . 

The program is based on the assumption that the hyperbolic arrival conditions at 
the planet are known. The hyperbolic excess velocity and the arrival asymptote for a 
given arrival Julian date are generally available (refs. 1 and 2) or can be generated with 
Lambert's theorem and the appropriate ephemerides. To establish an elliptical orbit 
about the planet, a deboost from a hyperbolic periapsis to an elliptical periapsis was 
assumed, and this assumption should not be restrictive for preliminary mission design. 
Occultation and imagery parameters within the orbit are then generated. The calculations 
have been greatly simplified by fixing the celestial bodies during one spacecraft revolu- 
tion. This simplification is valid since the angular displacement of the planet in its orbit 
about the Sun is small for a period of a few hours. A program option allows the investi- 
gator to examine future conditions by stepping the planet in time within its orbit about the 
Sun. The second zonal harmonic J 20 can be incorporated in these calculations to per- 
turb the argument of periapsis and the longitude of the ascending node of the satellite orbit 
by the appropriate input. 

In an attempt to make the program a useful analytical tool, emphasis has been 
placed on the speed and simplicity of the calculations . Keplerian mechanics and analytic 
solutions were used entirely throughout the program to make the calculations as rapid as 
possible. The mean heliocentric positions of Venus, Earth, and Mars are generated with 
mean orbital elements. The calculation sequence has been divided into a number of sepa- 
rate calculations which form the subroutines of the program. This modular form should 
facilitate the understanding of the many calculations performed. 

The necessary information pertaining to the implementation of the program is con- 
tained in the appendixes of this paper. The function of the main program is outlined with 
a flow diagram and the purpose of each subroutine is set forth (appendix A). A complete 
FORTRAN listing is included (appendix B) in addition to a sample input and output case 
(appendix C). 


SYMBOLS 

a semimajor axis, kilometers 

A,B angles defined by figure 5, degrees 
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AU astronomical unit, kilometers 

Co,Ci,C 2 ,C 3 ,C 4 coefficients of a polynomial 

d vector defined by figure 7 

d magnitude of d 

D Julian days since 1900 

D = — — — 

10 000 

e eccentricity 

E unit vector from planet toward a celestial body 

F,G,H,I angles defined in figure 2, degrees 

f true anomaly, degrees 

h a apoapsis altitude, kilometers 

hp periapsis altitude, kilometers 

i inclination, degrees 

i , j , k unit coordinate vectors 

JD Julian date, days 

J 20 second zonal harmonic of planet 

l semilatus rectum, kilometers 

M mean anomaly, degrees 

n mean angular rate, radians/second 

P,Q,R points in spherical triangle (see fig. 2) 
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P,Q,W unit coordinate vectors in PQW coordinate system; P points toward peri- 
apsis, Q is in orbital plane advanced to P by a right angle in direction 
of increasing true anomaly, and W completes right-handed system 

R rotational matrix from mean Earth equinox and equator of date coordinate 

system to mean planet equinox and equator of date coordinate system 

RPQW rotational matrix from XYZ TO PQW coordinate system 

r radius vector 

r magnitude of r 

r s surface radius of planet, kilometers 

S unit vector from center of planet parallel to arrival asymptote of incoming 

hyperbola 

5 magnitude of S 

t whole number of Julian years since 1950 

t' time past deboost, seconds 

T e Julian centuries since 1900 

V velocity, kilometers/second 

X,Y,Z rectangular Cartesian coordinates 

a right ascension of axis of rotation of planet referenced to mean Earth equinox 

and equator of date coordinate system, degrees 

j3 orbital plane orientation angle, degrees (see fig. 1) 

y declination of axis of rotation of planet referenced to mean Earth equinox and 

equator of data coordinate system, degrees 

6 declination, degrees 

A angle defined by figure 2, degrees 
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e mean obliquity of ecliptic, degrees 

A. right ascension, degrees 

li gravitational constant, kilometers^/ second^ 

r flOOTg], where [x] is fractional part of X, years 

T first point of Aries 

<f> angle between S and P, degrees 

P angle between planet-Sun vector and planet-satellite radius vector, degrees 

S2 longitude of ascending node, degrees 

co argument of periapsis, degrees 

Subscripts : 


C 

EC 

EE 

f 

h 

l 

P 

PE 


p,q,w 


px,py,pz 


denotes star Canopus 

measured in mean Earth equinox and ecliptic of date coordinate system 
measured in mean Earth equinox and equator of date coordinate system 
first value 

denotes hyperbolic orbit 
last value 

denotes planet Venus or Mars 

measured in mean planet equinox and equator of date coordinate system 

rectangular Cartesian components in PQW coordinate system 

rectangular Cartesian components in mean planet equinox and equator of date 
coordinate system 
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denotes S 

first and second root of an equation 
infinity 
denotes Sun 
denotes Venus 
denotes Earth 
denotes Mars 
initial conditions 

Dot over symbol denotes derivative with respect to time. 

METHOD OF CALCULATION 

The hyperbolic arrival conditions at the planet are inputs to the program. These 
inputs include such quantities as the arrival Julian date, the right ascension and declina- 
tion of the S-vector expressed in the Earth equatorial coordinate system and the hyper- 
bolic excess velocity. The S-vector is identical to the unit vector from the center of the 
planet parallel to the arrival asymptote of the incoming hyperbola. The size of the planet - 
ocentric ellipse must also be included in the form of the apoapsis and the periapsis alti- 
tudes. Since the orbital plane can be rotated about the S-vector with only small expendi- 
tures of fuel, the orientation of the orbital plane is considered to be an independent varia- 
ble and is at the discretion of the investigator. The plane is restricted, however, to con- 
tain the S-vector, but can assume any inclination that fulfills this requirement. The angle 
jS (fig. 1) has been taken as a measure of the orientation. The program allows for this 
angle to be incremented through any range of interest from 0° to 360°. For a given /3 
the computation of Keplerian orbital elements is based on a deboost from a hyperbolic 
periapsis to an elliptical periapsis. Imagery data and occultation data based on the posi- 
tions of the celestial bodies at the time of periapsis passage are then generated. 

Hyperbolic Arrival Conditions 

The hyperbolic arrival conditions at the planet for a given Julian date as obtained 
from references 1 and 2 are expressed in a coordinate system oriented to the mean Earth 


S 

1,2 


© 

? 

© 

o' 

0 
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Axis of rotation 



S (unit vector 
from center of 
planet parallel 
to approach 
asymptote of 
incoming hyperbola) 


Vector 
perpendicular to 
line of nodes in 
orbital plane and 
in direction of 
spacecraft motion 


Vector perpendicular to line 
of nodes directed easterly in 
equatorial plane 


Figure 1.- Orbital plane geometry, (p is measured counterclockwise in 
plane perpendicular to line of nodes.) 


equinox and equator. Since it is desirable to perform all the orbital calculations in a 
coordinate system oriented to the mean planet equinox and equator, the transformation 
between the two systems is necessary. 

Consider figure 2 where a. and y are the right ascension and declination of the 
planet's axis of rotation, respectively, expressed in the Earth coordinate system. For 
Venus (ref. 3, p. 20), 

a = 272°. 75 
7 $ = 71°. 50 

and for Mars (ref. 4, p. 334), 

= 317°. 793416667 + 0°.006520833t - 0°.001013 t 

= 54°. 6575000 + 0°.0035t - 0°.000631r 

where 

T = JD - 2 415 020 
e 36 525 

r = |100T^] where [X] is the fractional part of X 
t = 100T e - r - 50 
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North pole 
of Earth 



Figure 2.- Geometry of transformation from the Earth 
equator to the planet equator system. 


and JD is the Julian date at the time of inter 
est. The longitude of the ascending node of 
the Venus orbital plane on the ecliptic and the 
inclination of the Venus orbital plane to the 
ecliptic is (ref. 4, p. 113) 

£2^ = 75°46’46".73 + 3239” .46T e + l”.476T e 2 

i^ = 3°23'37”.07 + 3”.621T e - 0”.0035T e 2 

For Mars (ref. 4, p. 335) 

Go? = 48°47'll”.19 + 2775”. 57T e - 0”.005T e 2 

i^ = 1°51’01”.20 - 2”.430T e + 0”.0454T e 2 

The mean obliquity of the ecliptic e is given 
by (ref. 4, p. 98) 


e = 23°27'8".26 - 46".845T e - 0".0059T e 2 + 0".00181T e 3 


The rotation from the Earth equatorial system to the planet equatorial coordinate 
system is given in reference 4 (p. 330). It can be seen that this rotation is essentially 

an Euler rotation where the X-axis is rotated from T 
to Q along the Earth equatorial plane and followed by 
a rotation about Q through the angle 90° - y, and 
finally the X-axis is rotated along the planet's equato- 
rial plane to the ascending node of the orbital plane on 
its equatorial plane. The angle TQ is immediately 
seen to be a + 90°, and the rotation about Q is 
90° - y. The only angle remaining to be calculated is A which necessitates the solu- 
tion of two spherical triangles. Consider the triangle TQP as shown in figure 3. The 
angles F, G, and H are computed by 

cos H = cos £2 cos (90° + a) + sin £2 sin(90° + a)cos e 



Figure 3.- Spherical triangle defined 
in figure 2. 


sin H = \jl - cos 2 H 


cos F = 


cos (90° + oi)sin £2 - sin(90° + qi)cos £2 cos e 
sin H 


. sin e sin(90 + a) 
sin F = - 

sin H 
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ii i 


in it 


cos Q, sin(90° + a) - sin £2 cos(90° + a)cos e 


cos G = 


sin H 


sin G = sln ^^ 

sin H 

Now, consider the triangle PQR (fig. 4): The angle I is 
given by 

cos I = -cos(180° - F + i)cos(90° - G + y) 

+ sin(180° - F + i)sin(90° - G + y)cos H 
sin I = - cos2l 

and the angle A by 

sin(180° - F + i)sin H 

sin A = 

sin I 



Figure 4.- Angle relationships 
of spherical triangle defined 
in figure 2. 


COS A = 


cos(180° - F + i)sin(90° - G + y) + sin(180° - F + i)cos(90° - G + y)cos H 

sin I 


A = tan-i/^HLA^ 

\cos A J 

Therefore, the three Euler angles (90° + a, 90° - y, A + 180°) constitute the Euler rota- 
tion R which transforms the S from the Earth equatorial to the planet equatorial sys- 
tem, that is 

SpE = RS E e 

where 


r = 


cos (A + 180°)cos(90° + or) 

-cos(90° - y)sin(90° + a)sin(A + 180°) 

-sin(A + 180°)cos(90° + a) 

-cos(90° - y)sin(90° + a)cos(A + 180°) 

sin(90° - y)sin(90° + a) 


cos (A + 180°)sin(90° + a) 

+cos(90° - y)cos(90° + cn)sin(A + 180°) 

-sin(A + 180°)sin(90° + a) 

+cos(90° - y)cos(90° + a)cos(A + 180°) 

-sin(90° - y)cos(90° + a) 


sin(A + 180°)sin(90° - y) 
cos(A + 180°)sin(90° - y) 
cos (90° - y) 



Orbital Transfer to Planetocentric Ellipse and 
Orbital Perturbations With Time 

If the characteristics of the incoming hyperbolic orbit in the planet equatorial sys- 
tem (SpEpV^), the shape of the desired elliptical orbit about the planet (hp,h a j, and the 
orientation angle /3 are known, the Keplerian orbital elements can be determined for a 
transfer from a hyperbolic periapsis to an elliptical periapsis. If 

S PE = ( S px’ S py» S pz) 


then the right ascension and declination of the Sp E are, respectively. 

S r 


sin X 


5 py 


S,PE 




2 , c 2 
IX + S py 


cos X 


S,PE 


°px 

q 2 , c 2 

°px + °py 


X S,PE = tan 


-1 


/sin X 


S,PE 


and 


3 s,pe 


= sin 


cos X S;PE y 

‘ 1 ( S pz) 


(-90° 2 « s pE s 90°) 


The semimajor axis and the eccentricity of the elliptical orbit are given by 


a = 


2 r s + h a + hp 


r s + h a 1 
e = -2 1 


where r s is the surface radius of the planet. 

Consider figure 5. From figure 5 it can be seen that 

/tan 6 S pE \ 

A = sin - -M 1 

|tan y3 1 


(-90° ^ A ^ 90°) 


B = sin' 


/sin 6 s>pE ^ 


sin /3| 


(-90° ^ B ^ 90°) 
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The angle <£ between the vector S and the vec- 
tor to the periapsis P can be found by considering 
the geometry of the incoming hyperbola where the 
orbit is defined by 


r = ^ " e h 2 ) 
1 + e h cos f 


( a h <0; e h > 1) (1) 


As the radius r approaches infinity, the denomina- 
tor must approach zero, or 


1 + e h cos f^ = 0 

which can be solved for the cosine of the true anom- 
aly, 


From figure 6 it can be seen that 


z 



Figure 5.- Orbital transfer geometry. 


fco = 1800 - <#> 

so that 


cos f^ = -cos 0 


The two expressions for the cos f^ can be used to obtain 


cos d> = — 
V e h 


( 2 ) 


Evaluating equation (1) at a true anomaly of zero gives 

r s + h p = a h(! ~ e h) 


or 


e h = 1 


r s + h P 
a h 


( 3 ) 


The semimajor axis of the hyperbola can be expressed in terms 
of known parameters by evaluating the relation 



as r approaches infinity; that is, 



Figure 6.- Geometry of 
hyperbolic orbit. 


a h “ 


-M 


V 2 

v oo 


( 4 ) 
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Equation (2) can now be expressed in terms of known parameters by equations (3) and (4) 
as 

cos 0 = ^ (0° S(j) s 180°) 

p + (r g + hpjVoo 2 

The inclination defined between 0° and 180° at the ascending node is given by 

i = 0 (o° < js < 180°) 

i = 360° - p (l80° < p < 360°) 

The argument of periapsis and the longitude of the ascending node can be determined by 
considering the geometry; that is. 


or 


and 


or 


<jo = B - (p 


a ~ a s,pe " A 


(*) = B - 0 
£2 = Xg PE A 


w = 180° - B - cj) 

£2 = 180° + Xg ppj ~ A 

<jo = 180° - B - (p 
f2 = 180° + Xg pg + A 


(0° < p < 90°) 



< j8 < 180°) 


(l80° < p < 270°) 


(270° < p < 360°) 


Therefore, the planetocentric ellipse is defined by the Keplerian orbital elements a, e, 
i, (jo, and £2 for a transfer from one periapsis to the other. 

As the planet proceeds in its orbit about the Sun, the longitude of the ascending node 
and the argument of periapsis are perturbed by the second zonal harmonic J 2 q- These 
two orbital elements are expressed as functions of time as (ref. 5) 

£2 — £2 q + £2t 

Cl) = co Q + d)t' 
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where 



The elements £7 0 and oj q are the orbital elements at the time of deboost into the 
planetocentric ellipse and t' is the time past deboost in seconds. 

The gravitational constant and the radius of the planet are necessary in calculations 
involving the orbital elements and their perturbations with time. The gravitational con- 
stants for Venus and Mars are given by (ref. 6, p. 8) 

= 3.2485340 x 10 & km3/sec2 

= 4.297780 X 104 km3/sec2 
and the surface radius of Venus is (ref. 3, p. 18) 

r,, = 6085 km 

¥ 

and the surface radius of Mars is (ref. 7, p. 24) 

r S( y. = 3395 km 

Satellite Occupation of Celestial Bodies 

Before the positions in the orbit at which the satellite enters and exits the shadows 
of the celestial bodies can be computed, the unit vectors from the planet toward the celes- 
tial bodies must be found. These vectors are found by considering the mean orbital ele- 
ments of the planets about the Sun which lead to the heliocentric position vectors of the 
planets. The proper addition of these vectors will yield the desired vectors. Let 

D = JD - 2 415 020 


10 000 
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since the mean orbital elements are referenced to the mean equinox and ecliptic of 1900. 
The value of the astronomical unit in units of kilometers is given by (ref. 8, p. 6) 

AU = 149 598 845 km 

and the gravitational constant of the Sun is (ref. 6, p. 8) 

= 1.32715445 X 10 11 km 3 /sec 2 

Since the mean orbital elements of the Sun about the Earth (ref. 4, p. 98) are similar to 
the orbital elements of the Earth about the Sun, the semimajor axis is 

a ffi = 1. 00000023 AU 
© 

The eccentricity is 

e e = 0.01675104 - 0.00004180T e - 0.000000126T e 2 
The longitude of periapsis is 

u> 0 = 101°. 220833 + 0°.000047068D + 0°.0000339D 2 
and the mean anomaly is 

M e = 3 58°. 475845 + 0°.9856002670D - 0°.0000112D 2 - 0°.00000007D 3 

The inclination of the Earth's orbit is zero by definition, and the argument of the ascending 
node is taken as zero. The mean heliocentric position and velocity of the Earth are 
obtained by converting the orbital elements to Cartesian coordinates. The mean orbital 
elements for Venus are (ref. 4, p. 113) 

a^ = 0.7233316AU 

e<£ = 0.00682069 - 0.00004774T e + 0.000000091T e 2 

= 75°46’46".73 + 3239".46T e + l”.476T e 2 

= 130°9'49".8 + 5068". 93T e - 3".515T e 2 - ^ 

i^ = 3°23 , 37".07 + 3".621T e - 0".0035T e 2 

M<j) = 212°.603219 + 1°.602130 1 540D + 0°.000096400D 2 

and for Mars (ref. 4, p. 113) 

SLtf = 1. 523691 5AU 

e^ = 0.09331290 + 0.000092064T e - 0.000000077T e 2 

Slrf = 48°47' ll".19 + 2775". 57T e - 0".005T e 2 - 0".0192T e 3 

c Oq, = 334°13'05".53 + 6626". 73T e + 0".4675T e 2 - 0".0043T e 3 - 


I 
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= l°5l'l".20 - 2”.430T e + 0".0454T e 2 

= 319°. 529425 + 0°.5240207666D + 0°.000013553D 2 + 0°. 00000002 5D 3 

Therefore, the mean heliocentric, ecliptic positions of the Earth and the planet (that is, 
F 0e and r Q p) can be obtained. 

The vector from the planet to the Earth is 

F P® = r o® - ?®p 

and the vector from the planet to the Sun is 

F Po = - F oP 

These vectors are rotated from the ecliptic to the Earth equatorial system by 


X EE 


1 0 0 


X EC 

y ee 

= 

0 cos e -sin e 


y ec 

_ Z EE 


0 sin 6 cos e 


Z EC 


and to the planet equatorial system by the Euler rotation previously defined; that is, 


X PE 



X EE 

y pe 

= 

R 

y ee 

Z PE 



Z EE 


The unit vectors to the celestial bodies are therefore 

E = (E px ,E py ,E pz ) 

where 

r PE = ( X PE 2 + y PE 2 + Z PE 2 ) ^ 

The unit vector to the star Canopus in the planet equatorial system is found by first 
rotating the vector from the mean Earth equinox and equator of 1950 coordinate system 

X c = -0.060340592 
Y c = 0.60342839 
Z c = -0.79513092 
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to the mean Earth equinox and equator of date. This rotation is described in reference 4 
(pp. 28 to 39). The updated position vector is then rotated to the planet equatorial system 
by 

X C,PE = RX C,EE 

The solution to the occultation problem is from reference 8 (p. 155) and the sim- 
plified geometry is shown in figure 7. The unit vector toward the body being occulted is 
defined by 

E = E px f + E py f + E pz k 


which can be rotated to the PQW coordinate system; that is, 

E = EpP + EqQ + E w W 

The radius vector to the spacecraft r" can be expressed in the PQW system by 

r = r cos fP + r sin fQ + OW 


where 

a(l - e 2 ) i 

I* — — 1 : — 

1 + e cos f 1 + e cos f 


(5) 

( 6 ) 

(7) 


The geometric constraint can be obtained from the condition that upon entrance to or exit 
from the shadow, d is parallel to E; therefore, 

E • d = -d 


where 


and 


d = r - r s 
d = \/r 2 - r s 2 



Figure 7.- Geometry of occupation. 


Therefore, 

E • (r - r s ) = -\p - r s 2 
and since E is perpendicular to Fg, 

E • r = ~\p- - r s 2 (8) 

Expanding the dot product with equations (5) 
and (6) and eliminating r with equation (7) 
allows equation (8) to be expressed in terms 
of only the true anomaly; that is, 
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( 9 ) 


Ep£ cos f EqZ sin f 
1 + e cos f + 1 + e cos f 


£ 2 

(l + e cos f)^ 



1/2 


which can be reduced to standard form as 

Cq cos4f + Ci cos3f + C2 cos2f + C3 cos f + C4 = 0 


where 

c ° 'ft )" 64 - 2 (t) 2 (V ‘ E p 2 ) e 2 + (V + V ) 2 

Cl - - 4 (t) 2 (V - E P 2 ) 6 

C 2 . .^)\* - 2 (^) 2 ( Eq 2 . Ep 2) - $)* (' - V> 2 

+ 2(E q 2 - E p 2 )(l - E q 2 ) - 4E p 2 E q 2 
c 3 - 4(^) 4 e - 4(^) 2 (l - E q 2 )e 


All the real roots of equation (9) are extracted by Descartes technique (ref. 8, p. 430). 
The spurious roots are then rejected by 



= r 


s 


r • E < 0 


It can be shown that the spacecraft enters the shadow at fi provided that 

2r s 2 (l + e cos fi)(-e sin fi) + 2Z 2 (E p cos fi + E q sin fi)(-E p sin fi + E q cos f^> 0 

where fi is a solution of equation (9). The exit from the shadow is at f 2 . Once fi 
and f 2 have been determined, pertinent parameters such as the length of time in the 
shadow can be computed. 


Vertical Photography at a Specific Lighting Angle 

The positions in the orbit of the satellite about the planet that correspond to desira- 
ble lighting conditions on the surface for vertical photography are of interest in mission 
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analysis. These positions along with the orbital elements lead to the definition of the sur- 
face which can be photographed. 

The lighting angle \p is the angle between' the planet-Sun vector and the planet- 
satellite radius vector. When the lighting angle is 90°, the satellite is directly over the 
terminator. The unit vector from the planet toward the Sun is rotated to the PQW coor- 
dinate system by 


Ep 




E px 

E q 

= 

RPQW 


E py 

E w 




E pz 


where RPQW is the rotation matrix from the planet equatorial to the PQW system and 
is defined by 

cos co sin £2 
+sin go cos i cos £2 


RPQW = 


cos co cos £2 
-sin co cos i sin £2 


sm co sin l 


-sin co cos £2 
-cos co cos i sin £2 

sin £2 sin i 


-sin co sin £2 
+cos £2 cos co cos i 

-sin i cos £2 


cos co sin i 


cos 1 


The unit vector from the planet toward the spacecraft in the PQW coordinate system is 
given by 

r = cos fP + sin fQ + 0W 

Since ip is the angle between the Sun vector and the radius vector 

E • r = cos ip 


From the dot product 


or 


I E x r | = sin ip 
Ep cos f + Eq sin f = cos ip 
cos \p - Eq sin f 


cos f = 


E r 


From the cross product 

E w 2 sin 2 f + E w 2 cos 2 f + (E q cos f - E p sin f j 2 = sin2 ip 


( 10 ) 


(ID 
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Using equation (10) to eliminate the cos f in equation (11) and putting in standard form 
gives 

Cq sin f + Cj sin f + C 2 = 0 (12) 


where 


c 0 



2E q 2 + E p 2 


-2E a 3 

Cl = ^ - cos xp - 2Eq cos ip 

E p 2 


E — 

c 2 = E w 2 + cos 2 i// - sin 2 i// 


E r 


The quadratic formula is used to obtain the sines of fj and f 2 , and equation (10) is used 
to obtain the corresponding cosines. Therefore, 


f l,2 = 


tan" 


1 /sin i\2 \ 

\ cos *1,2/ 


There are two true anomalies that correspond to a given lighting condition. It is very 
important to determine whether the satellite is entering the desirable lighting conditions 
or leaving these conditions. Differentiating equation (10) with respect to time yields 


and since 


f > 0, 


!// = f 1 


(Eq cos f - Ep sin f 
^ -sin i// 


sign ( 1 //) = sign 


(Eq cos f - Ep sin f 


-sin 1 // 


If sign(i//) is positive, the lighting angle is increasing; if signO//) is negative, the 
lighting angle is decreasing. Distinction between ascending and descending motion at the 
true anomalies in question is made by the following: 

Ascending motion: 

-90° < u> + f < 90° 


Descending motion: 


90° < a> + f < 270° 
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RESULTS AND DISCUSSION 


The computer program has three distinct advantages. The use of Keplerian mechan- 
ics and analytic solutions to the occultation and imagery problems has resulted in an 
extremely fast program. The use of mean orbital elements to generate the ephemerides 
of the planets is considered to be an advantage in that the program is self-contained and 
does not require the use of ephemeris tapes. Possibly, the major advantage of the pro- 
gram is the modular form in which it is written. This form allows the user to insert 
additional calculations into the program with a minimum of effort. The subroutines 
(appendix B) have been written in as general a form as possible and can be easily adapted 
to various orbital- mechanics problems. The modular form in which VAMOOS was 
written resulted in a field length (storage requirement) of 40 000g. 

Sample plots which have been found useful in design studies for a Mars mission are 
presented in figures 8, 9, and 10. The orbital inclinations which result in Sun, Earth, or 
Canopus occultations during the first revolution of the satellite about the planet can be 
seen in figure 8. In addition, the duration of the occultation period is given. This type 
of plot has been found useful in selecting an orbital inclination that insures visibility of 
the satellite from Earth and provides adequate sunlight for the satellite. The computing 



Figure 8.- Sun, Earth, and Canopus occultation durations during first satellite revolution 

about Mars. 
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time on the digital computer for the data presented in figure 8 was 0.16 second. Ground- 
track and lighting-angle plots (fig. 9) have been found useful in examining the imagery 
requirements for an orbital planetary mission. For a given orbital inclination the 
sequence of events during many satellite revolutions can be generated with the program. 
(See fig. 10.) 


CONCLUDING REMARKS 

A program for the calculation of mission design data for either Mars or Venus has 
been developed. The planetocentric elliptical orbit about the planet is examined to deter- 
mine the times at which the satellite enters and exits the shadows of the Sun, Earth, and 
Canopus. The photographic coverage of the surface of the planet is also determined. 
Because of the speed and simplicity of the calculations, VAMOOS has proved to be a use- 
ful and flexible tool for mission design. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., December 13, 1968, 

194-82-01-05-23. 
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PROGRAM DESCRIPTION 


The VAMOOS program has been written entirely in FORTRAN IV computer language 
for the Control Data 6600 digital computer and contains a main program and 27 subrou- 
tines. The modular form in which the program was written resulted in a field length 
(storage requirement) of 40 OOOg. Very few calculations are performed in the main pro- 
gram since its primary purpose is one of administration. The input data are read into 
the main program and the appropriate subroutines are used to perform the desired com- 
putations. The angle /3 (fig. 1) can be incremented from 0° to 360° in whole degrees at 
the discretion of the user. An additional program option allows for the range of j3 to be 
investigated at a future date by stepping the planet in time within its orbit about the Sun. 
This option permits the time after arrival to be incremented over any time period with 
printout as frequently as desired. The second zonal harmonic J 20 can be incorporated 
in these calculations to perturb u> and 12 of the satellite orbit by the appropriate input 
if such is desired. All these options are controlled by the main program which is out- 
lined in the flow diagram. The letters in parentheses correspond to the subroutine that 
performs the stated calculation. A brief statement of the purpose of each subroutine con- 
tained in the program follows : 


JULCAL 


Converts Julian data to calendar data 


REQVEQ 

REQMEQ 

REQPEQ 

RECEQ 

RXYZPQW 

RPQWXYZ 


Rotates a vector from the mean Earth equinox and equator of date 
coordinate system to the mean Venus equinox and equator of date 
coordinate system 

Rotates a vector from the mean Earth equinox and equator of date 
coordinate system to the mean Mars equinox and equator of date 
coordinate system 

Rotates a vector from the mean Earth equinox and equator of date 
coordinate system to the mean planet equinox and equator of date 
coordinate system 

Rotates a vector from the mean equinox and ecliptic of date coordi- 
nate system to the mean Earth equinox and equator of date coordi- 
nate system 

Rotates a vector from the XYZ to the PQW coordinate system 

Rotates a vector from the PQW to the XYZ coordinate system 
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E EARTH 
EVENUS 
EMARS 
SUN 

SUNBAND 

OCCULT 

CONCAR 

CARSPH 

TCONIC 

TINVS 

QADRAT 

CUBIC 

QARTIC 

CROSS 

DOT 

ORBIT 

PRECES 

LATLNG 


Calculates the mean heliocentric position and velocity of Earth 

Calculates the mean heliocentric position and velocity of Venus 

Calculates the mean heliocentric position and velocity of Mars 

Calculates the positions in orbit which correspond to various lighting 
angles and writes output data 

Calculates the two positions in orbit which correspond to a given 
lighting angle 

Calculates the entrance and exit true anomalies of occultation 
Converts conic elements to Cartesian coordinates 
Converts Cartesian to spherical coordinates 

Calculates the time from periapsis passage for a given true anomaly 

Converts mean anomaly to eccentric and true anomaly 

Solves the equation AX^ + BX + C = 0 for the real roots 

Solves the equation AX 3 + BX2 + CX + D = 0 for the real roots 

Solves the equation AX^ + BX3 + CX2 + DX + E = 0 for the real roots 

Calculates the vector cross product 

Calculates the angle between two vectors 

Calculates Keplerian orbital elements for a periapsis to periapsis 
deboost 

Transforms mean Earth equinox and equator coordinates from one 
epoch to another epoch 

Converts a Cartesian position to latitude and longitude 
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EULER 

VECTOR 


Performs an Euler rotation 

Calculates the position of the Sun, Earth, and Canopus in the planeto- 
centric, planet equator, coordinate system and writes output data 
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planet equator (REQVEQ or REQMEQ) 



4 

Perturb orbital elements u> and 
S2 by J 20 if TIME > 0 
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PROGRAM LISTING 


A complete FORTRAN IV listing of the program VAMOOS is contained herein. The 
purpose of each subroutine is stated in the listing along with a brief definition of the input 
and output parameters. If additional subroutines are required for the computations con- 
tained in a specific subroutine, this fact is also designated. 


C 

C VAMCCS PROGRAM - VENUS AND PARS ORBITAL OCCILTATION SIMULATOR 

C PRO GRAMM t D FOR THE CDC 6600 

C 

C THIS PROGRAM CALLS SUBROUTINES JULCAL, REQVEU, REQMEO, REOPEO, RECEQ, 

C RXYZPCW, RPCWXYZ, EEARTh, EVENUS, EMARS, SUN, SUNBAND, OCCULT, CONCAR, 

C CAR SPF , TCCNIC, TIN\S, CACRAT, CUBIC, CART I C , CROSS, DOT, ORBIT, 

C PRECES, LATLNG, EULER, VECTOR 

C 

C 

DIMENSION G ( 6 ) ,E(6) ,PPGW ( 3 , 3 ) , SUNB ( 6 ) 

REAL JD 

NAMELIST /CASE/ JC , SVCEC E , S VRAE , VINF , H A , HP , BF F S T , B L A ST , B STE P , T L A ST , T 
1 STEP , X J2C , I EOCY ,SLN1 ,SUN2 ,SUN3 , SUN4 , SUN5 , SUN6 
ANGLE! X >=AMCD!X, 3 60.) + 180. -SIGN! 1 80 . , X ) 

( 

RD=57.295779513C£23 
DR =C .0174532^251994 
PI=3. 14159265358979 
UVENUS=324653.4 
R V ENL S=6C 85 . 

UMARS=42977.8 

RMARS=3395. 

C 

6 DC TL AS T =0 • 

T S TL P=C • 

BFR S T =C . 

B S 1 E P =C • 

B L A S T = 0 . 

T I ME =C . 

X J 20 =0 . 

SVUECE=0. 

S V R A E = 0 . 

SUN 1 = 0 . 

SU N2 =0 • 

SUN 3 = 0. 

SLN4=C • 

SUN5 = 0 • 

SUN 6=0 . 

C 

R£AC(5,CASE) 

WRITE (6, CASE) 

I F ( ECF , 6 ) 8C0 * 1 
C 
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1 CCNTINUE 
SUNB ( 1)=SUN1 
SLNB ( 2) =SUN 2 
SLNB ( 3 ) = SUN3 
SUNB! A )=SUN A 
SLNE ! b ) = SUN5 
SUNB! 6 )= SUN6 

( 

IF ( I BODY. EO. 2) GO TO 20 
IF! IE0DY..EQ.A) GO TC 21 
C 

20 U=LVENLS 
R S=R V ENL S 
GO TO 22 

21 U = IMRS 
RS=RMARS 

C 

22 CONTINUE 

W J C= FLO AT ( I F I X ( JO) ) 

F JD= JD-WJD 

IF ( F JD .EG • . 5 ) F JL)=F JD«-0.0000 1 
CALL JULC AL ! G *W JD, F JD, C ) 

I YA= IF I X ( G ( 1 ) ) 

IMA= IF I X ! G l 2 ) ) 

ILA= IF I X( G C 3 ) ) 

C 

2 CALL JLLC AL!B,WJD,F JD,0) 

IYB= I F I X ( B< 1 ) I 

IRB= IFI X(B (2) ) 

ICB= 1 F I X ( B ( 3 ) ) 

C 

WRITE (6, 101) JD.IYE ,IM E , I C E , IY A, IMA, IDA 
C 

CALL VECTOR! JC.C ECS ,R AS , DEC E,RAE, DECC, RAC, S>, SY.SZ, EX, EY,EZ ,CX,CY, 
1CZ , IB JD Y ) 

C 

IF! T I ME .GT. .0 1) GO TO 7 
IF! IBODY .EG .2 ) GC TO 3 
IF < IBCLY.EU.A ) GC TC A 
C 

3 CALL REGVEC < J C , COS ( SVD ECE*CR )*COS(SVRAE*DR) ,COS ( S V DEC E*DR ) *SI N l S VR 
1AE*LR)«SIN( SVDEC E #D P ) , X SVP ,Y S VP , Z S VP , S VDEC P ,SVRAP) 

GO TO 5 

A CALL RECMEC!JD,CCS!SVDECE*CR )*COS( S VRAE *DR ) ,COS! SVDEC E*DR ) *S I N ( SVR 
1AE*GR ) ,SIM SVUECE*D R) , X SVP ,YS VP , Z SV P , SVDEC P ,SVRAP ) 

C 

5 IFIBSTEP.LT..001) GC TC 6 
1B1= IFIXt BFRST ) 

I B2 = I F I X ( BLAST) 

IB3= IF IXIBSTEP) 

GO TO 7 

6 181=1000 
IB 2= 1 CC 0 
IB3=1CCC 

C 

7 DC 11 1=1 El ,IE2,IE3 
C 
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BET A=PLOAT (I) 

IP ( I .EC. I COC ) BE T A= BFR ST 

IF (A6S(SIN(BETA*QR) ) .LT . AB S I SI N I SVDECP *DR ) ) ) GO TC 11 
C 

CALL CRB I T< SVC ECP , S \.PAP ,V INF ,H A, HP, BET A, A , E ,X I , W , 0 , POEC ,PRA ,U,RS ) 
IFITIME.LT. .011 GC TO 8 
XN=SCRT (U/A**3) 

0D0T=-3 ./2 . *XJ20*XN* (RS/ A ) ** 2/ ( 1 ,-E*E )** 2 *CCS( XI *DR ) 

WD0T=3. *X J2C*XN*< RS/ A) **2/(l.-E*E)**2*(l.-5 ./4 .*S INIX I *DR) ** 2 ) 

W = W + I MiC 0 T * T IME+24.*3£OC . ) *RD 
G=C+( GDCT*TIME*24.*3600. >*RC 
PDEC = AS IN< SINIW*DR)*SIN (XI*DR) ) *RD 
SSDEL=T AN(FCEC*CR)/TAN (X I*CR ) 

CSDEL=SSDEL/TAN(W*DR)/COS(XI*DR) 

SOEL = AT AN 21 SSDEL*CSDEL ) *RD 
PR A= ANGLE (O+SOEL) 

8 CONT INUE 
C 

IF ( IBC0Y.EQ.2 ) WRITE(6,102) BETA, TIME 

IFI IBOCY.EQ .4 ) WRITEI6 * 1C 3 ) BETA, TIME 
C 

CALL RXYZPCW (0. ,C. ,0. »X I ,Vn ,0 ,RPOW ,VP ,VQ, VW ) 

CALL COT ( +RPQWI 2, 1 ) , -*RPCW (2,2) ,+RPUW(2 ,3) , S X , SY , SZ , THRUST ) 

CALL DLT I COS I SVDECP *CR ) *CCS I SVRAP*DK ) CO S ( SVDECP*DR)*SIN(SVRAP*UR) 
1 ,SINI SVDECP*DR) , SX ,SY , S Z , Z A P ) 

SVP = AC0SIU/(U-MRS-*FP)*VINF**2))*R0 


V PE = SORT I U*( R S + P A )/A/(RS + HP) ) 

VPH=Si)RT(VINF**2 + 2.*U/ ( PS + HP)) 

DELV=VPH-VPE 

VA l-=VPE*( (1 .-E ) / C 1 . ■* E ) ) 

PEPIOD=2.*PI*SORT ( A**3/U)/36C0. 

C 

V>R IT E C 6 ,10 A ) SVDECE , SV P *E , V I N F ,S VDE C P , S VR AP ,H A, HP , J D , PDEC , P P A , ( (RP 
lQkIJ.K) ,K = 1 , 3 ) , J = 1 , 3 ) .THRLST , A , E , X I , W , 0 , DEC S , R AS , OECE , R AE , SX ,SY,SZ 
2»E>»EY»FZ»PERICD*CX ,CY ,CZ,ZAP ,S V P , V PE , V AE , V FH , DEL V ,X J 20 
C 

IF ISUN1.LT. .0 1) GO TO 12 

CALL SUN! SX,SY,SZ,SLNB,A,E,XI,W,Q,RPQW,U,PS,1) 

C 

12 CALL CCCULTI A,E ,XI , V-,0,L,RS ,SX ,SY ,SZ,TOCC,T 1, ALT1, FI, DEC1, R A1,T2, A 
1LT2»F2,DEC2»RA2 ,KK ) 

WRITEIfc ,1C5 ) TCCC.Tl ,F1 , ALT 1 , DEC 1 ,R A 1 , T 2, F 2 ,ALT2, DEC2, RA2 
CALL OCCULT (A ,E ,XI ,W , □ , U, RS , E X , EY , E Z , TGCC ,T 1 , ALT 1, FI ,DEC1 , P A1 , T2 » A 
1LT2»F2»DEC2»RA2,KK) 

WRITE (6 ,1C6) TUCC ,T1,F1 ,ALT1 ,DEC1 , R A1 , T 2 , F 2 , ALT2 , DEC2 , RA2 
CALL OCCULT I A, E ,X I , W ,0 , U , RS ,C X , C Y ,C Z , TOCC ,T 1 , ALT 1 , F 1 , DEC 1 , RA1 , T2 , A 
1LT2 »F2 »DEC2,RA2 ,KK) 

WRITEU.1C7) TOCC.T 1 ,F 1 , ALT 1 ,DEC 1 , R A 1 , T2, F2 , ALT2 , DEC2 ,RA2 
C 

11 CONT INUE 
t 

IFITSTEP.LT..OC 1) GO TO 6CC 
IFIT 1ME.GE.TLAST) GC TO 600 
TIME=TIME+TSTEP 
WJC=toJO+FLOAT II FIX ITST EP) ) 

F JC=F JD + AMGDITSTEP ,1.) 

JC=W JC+FJD 
GC TC 2 


31 



n a o 


APPENDIX B 


c 

101 FORM £ T ( 1H l ////17H JULIAN DATE F2C. 1//22X ,22HYt A R MC 

1NTF UAY//22H CALENCAR CATE , I A , 5 X, I 3 ,7X , I 3 //22H ARRIVAL 

2 CAT E , I4.5X.I3 ,7X,I 3) 

102 FORMAT! #1 *////14X »* BETA =* » F9 .4 ,4X ,* ( VENUS EQUATOR, 

1VENUS E0UINCX)*/14X ,* TIME PAST ARRIVAL = *,F8.1,* DAYS*) 

IC3 F0FMAT(*1*////14X,* BETA =*, F9 .4, 4X ,*( MARS EQUATOR, M 

1ARS ECU INOX )* / 14X ,* TIME PAST ARRIVAL = *,F8.1,* DAYS*) 

104 FOFMAT (1HC/.9F SVC ECE=E 16 .8 , 10H SVRAE =E16.8,1CH VINF =E16.8 
ltlQH SVCECP=E16.8/9H SVPAP =E16.8,19H FA =E16.8,10H HP 

2 = E 1 6 . 8. 1CH JO =E 1 6 . 8 / SH PLAT =E16.8,10H P LONG_ = E16 . 8, 1 0 

3H PX = E 16 . 8 , 1 OH PY =E16.8/9H PZ = E16.8,10H CX 

4E16.8.10H OY =E16.8,10H QZ =E16.8/9H WX =E16.8,10H 

5 WV =E 1 6 • 8 » 1 CH WZ =E16.8,10H VEL SLN=E16 .8/9H SMA =E16 

6. 8, IGF ECC =E 16 .6, 10H INC = E16.8,10F ARGPER=E16. 8/9H AR 

7GNC0=EI6.8» 10 H CECSJN=E 16. 8, 10H RA SUN = E16.8,10H DEC ET H=E1 6 • 

88/9E RAETH =E16.8,10H XSUN =E16.8,1CH YSUN =E16.8,1CH ZSU 
9N =E16.8/SF X E ART F=E 1 6 . 8 , 10H Y E ART H = E16 .8 , 10H ZEARTH=E16 .8 , 1 
1 CH PERI0D=E16.8/SF XCANPS=E 16. 8 , 10H YC ANPS=E 16. 8 , 1 OH ZCANPS 

2=E16.£,1CH ZAP =E16.8/9H SVP =E15.8,10H VPE =E16.8,10H 

3 VAE =E16.8,1QH VPH =E16.8/9H DEL V =E16.8,10H XJ20 =E1 

46.8) 

1C5 F0RMATI22HC SLN OCCLLTATICN T I ME 1 5X , F 1 1 . 2 /2 3H ENTER SUN OCCULTATI 
10N.14X.5F11.2/2-2H EXIT SUN OCCUL T AT ION , 1 5X , 5F 1 1 .2 ) 

1C 6 FOFMAT (24HG EARTF OCCULTATIGN T I ME , 1 3X , F 1 1 . 2/2 5H ENTER EARTH OCCU 
1LT AT ION ,12X, 5F1 1 .2/24H EXIT EARTH OCCULT AT ION, 1 3X , 5F1 1 . 2 I 
1C 7 FORMAT (26HC CANOPUS OCCULTATION T I ME ,1 1 X ,F1 1. 2/27H ENTER CANOPUS 
10CCULTATI0N ,10X , 5F 1 1.2/26H EXIT CANOPUS CC CULTATI ON , 1 IX ,5 FI 1 . 2 ) 

C 

800 STOP 
END 

SUBROUTINE JULC AL( X ,W0 I , FD I, IND ) 

C 

C THIS SUBROUTINE CONVERTS A GIVEN JULIAN DATE OR THE NUMBER OF 

C WHOLE AND FRACTIONAL DAYS SINCE JANUARY 1, 1950, 0 HRS., TC THE 

C CO FF ESFCNC ING CALENEAR OATE. 

WDI - INTEGRAL PART OF JULIAN DATE GR WHOLE NUMBER OF DAYS SINCE 
JANUARY 1, 1550, C FRS. 

C FCI - FR ACT 1CNA L PART OF JULIAN DATE OR FRACTIONAL NUMBER OF DAYS 

C SINCE JANUARY 1, 1950, 0 HRS. 

C IND - CCNTRCL INTEGER. 0 IMPLIES JULIAN DATE, 1 INPLIES DAYS 

C XI1-6) - CALENDAR CATE ( Y EAR , MONTH , D AY , HOUR , MI NUTE , SECOND ) 

DIMENSION X (6 ) » A ( 12 > »W ( 12 ) 

WD = WD I 
FD=FDI 

I FI IND) 1,1 ,5 

1 IF< Ft- .5 >2, 2, 3 

2 FD=FU+. 5 
WD =WD- 1 . 

GO 10 4 

3 FD =FD- . 5 

4 WD=WD-24332 82 . 

5 WD=VD+1 • 

DY-365. 

Z = 2 . 

N = C 
Q=4 • 

6 WD=WD-DY 
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IF iWC) 10, 10,7 

7 N=N + I 
Z=Z*1. 

CK=0-Z 
1F(CK >9,9,8 

8 DY=3 65. 

FC = 2 8 • 

GO TC 6 

S DY = 26 6 . 

Q=C+A. 

FC =2 S. 

GO TO 6 

10 WC=VD+DY 

CO 11 1=1,12 

11 A ( I > =C • 

C 1 =3 1 . 

C2= 2 C * 

DO 13 1=1,12 
A ( 1 > = 1 • 

CA=FC*A(2>+C1*(4( 1 ) 4 A ( 3)+A(5)+AI7)+A<8)+A(lC)+A( 12)I+C2*(A(4>+ 
1 M6 )+AI9 1 + M 111) 

Vil I ) = rtD-C A 
IF(M I) >12,12,13 

12 IFU-1J 15,15,16 

15 M0N= 1 

GO TO 1 A 

16 MGN= I -1 
WD = fc( MC N ) 

MC N = MO N +1 
GO 70 14 

13 CONT INUE 

14 N = N + 50 
xt i>=n 

X ( 2 ) = NCN 
X ( 2) = fc\D 
FH=FC*24. 

N = F F 
X< 4 >=N 

FM = ( FH-X < A ) )*60. 

N = FK 
X ( 5 ) =N 

X ( 6 ) = l FM— X ( 5 ) ) *60 • 

RETURN 

END 


SI E ROUT IN E PEQV EQ t JD , XE C , Y E C , ZEQ , XV EQ , Y VEQ , Z VEQ , DEC VEQ , R A V E C ) 

C 

C THIS SUBROUTINE ROTATES A VECTOR FRCM THE MEAN EARTH EQUATOR- 

C EQUINOX TO THE MEAN VENUS E QU AT OR- ECU I NGX CCORDINATE SYSTEM. THIS 

C ROUTINE CALLS SUBROUTINES REQPEQ AND LATLNG. 

C 

t JD - JULIAN DATE AT TIME OF INTEREST. 

C XE C , Y EC , Z EQ - VECTOR IN THE EARTH EQUATORIAL SYSTEM 

C XVE C , YVEQ , Z VEQ - VECTOR IN THE VENUS EQUATORIAL SYSTEM 

C CECVEQ,RAVEQ - DECLINATION AND RIGHT ASCENSION OF THE VECTOR IN 

C THE VENUS ECUAT0RI4L SYSTEM 

C 

REAL JD 
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c 

TE = ( JO-24 15C20. ) /36f 25. 

C 

AL FH A0= 27 2.7 5 
GAFM4C=71 .50 

OMEGA=75.779f 4E+C . 699850* TE+0. 00041 l*TE**2 
>1=3 .3536 30+0.00 IOC 5*T E-C .5 7CE-6*TE **2 
C 

CALL RE CP EC ( JC .ALPHA C, GAM MAO. CM EGA* XI* XEQ, Y EC , ZE Q , XVE Q , YVEG , Z VEQ ) 
CALL LATLNG (XVEG»YVEG»ZVEC»DECVEQ»RAVEQ) 

C 

REURN 

ENC 


SUBRCUT I NE PECMEC ( JC ,X EG ,Y EQ*Z EG ,XMEO, YME C , ZMEG, DECM EQ, R AM EC ) 

C 

C THIS SUBROUTINE ROTATES A VECTOR FRCM THE MEAN EARTH EQUATCP- 

C EQUINOX TO THE MEAN MARS EQUATOR-EQUINOX COORDINATE SYSTEM. THIS 

C ROUTINE CALLS SUBROUTINES REQPEQ ANC LATLNG. 

C 

C JD - JULIAN CATE- AT TIME OF INTEREST 

C XEC.YEC.ZEO - VECTOR IN THE EARTH EQUATORIAL SYSTEM 

C XMEC .YMEQ ,ZMEO - VECTOR IN THE MARS EQUATORIAL SYSTEM 

C DECMEC, RAMEC - CECLINATICN AND RIGHT ASCENSION OF THE VECTOR IN 

C THE MARS EQUATOR I AL SYSTEM 

C 

PEAL JD 
C 

T E = I JC - 24 15C20. )/3£ 525. 

TAU=AMC0(TE*1C0.,1.) 

TP = 7E* ICC .-TAU-5C. 

C 

ALPFAC=3I 7. 79341 6667+0. 65 2C833E-2*TF-0.nriC 13 *TAU 
GAMMA0=54 .t5T5CCGCC+C.CC35*TP-0.000£31*TAU 
CM ECA=46.76644167+C. 77099167*1 E- 0.1 3886 889E -5 *TE**2 
X I = 1 .E5C333333-C.675E-3*TE+0.12611111E-4*TE**2 
C 

CALL RE CPECI JC , AL PF AO , C AMM AO, CMEGA * XI , XEC,Y EC , ZEO ,X MEQ, YME C, ZMEQ ) 
CALL LATLNG (XMEC , YM EC , ZMEQ , DEC MEQ, RAMEOI 
C 

PE TURN 
ENL 


SUBPCUT 1NE PECPEQ(JC,ALPHAO,GAMMAO,CMEGA, XI ,XEQ, YEQ,ZEQ,XPEQ , YPEQ , 
1ZPEQ) 

C 

TFIS SUBROUTINE ROTATES A VECTOR FRCM MEAN EARTH EQUATOR-EQUINOX 
TO PLANET EOUAT GR-ECUINCX COORDINATE SYSTEM. THIS ROUTINE CALLS 
SLBP CUT INE EULER. 
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t JC - J'JLIAN CATE AT TIME OF INTEREST 

C ALPFAO.GAMMAO - RIGHT ASCENSION AND DECLINATION OF THE PLANETS 

C AXIS CF ROTATION EXPRESSED IN THE EARTH EQUATORIAL 

C COORDINATE SYSTEM 

C OMEGA, XI - LONGITDCE 3F THE ASCENDING NODE AND INCLINATION OF THE 

C FLANETS ORBITAL PLANE REFERENCED TO THE ECLIPTIC AND 

C VERNAL EQUINOX 

C XEC.YEC.ZEC - COMPONENTS CF THE VECTOR IN TFE EARTH EQUATORIAL 

C COORDINATE SYSTEM 

C XPEQ.YPEO.ZPEQ - COMPONENTS OF THE VECTOR IN THE PLANET EQUATORIAL 

C COORDINATE SYSTEM 

REAL JC 

DIMENSION PPQW(3,3> 

DR=0. Cl 74 5 3 292 5 1994 3 
RC=57.255775513C823 
TE = ( JC-241502C. >/36 5 25. 

E = 23 .4 522 9444-0. 13C 125 E- 1*T E-0 . 1 63 6 8 88 9E -5 * IE **2+ 0. 50 27777 8 E- 6*TE* 

1*3 

C 

CE=COS( E*DR) 

SE = S INI E*DR I 
CAL=COS( ALPHAO*CR) 

SAL=SIN( ALPHAC*DR) 

CGM = CDS ( GAMMAG*DR I 
SGM=SIN(GAMMAC*CR) 

CCM=CCS(CMEG A*CR) 

SOM=SIN(OMEGA*DR) 

C 

CZP=CE*SCR*CAL-CCM*S AL 
SZ P = SQRT ( 1 .-CZP*CZF ) 

ZP=ATAN2ISZP,CZP)*RC 

SXP=SE*CAL/SZP 

CXP=(-CE*CCR*CAL-SCK*SAL)/SZP 

XP=ATAN2ISXP,CXP)*RC 

SYF-=SE*SUM/SZP 

CYF= ICE*SCM*S AL+CCM*CAL I/SZP 
YP=ATAN2( SYF,CYP)*RC 
C 

CI=CCS<<XP-XI)*CP)*SINIIY F-GAMMAC )*DR )+SIN( (XP-XI) *DR ) *COS ( I YP-GAM 
1M AO ) *DR )*CZ P 
SI=SCRT (l.-CI*C I ) 

C 

SNP=SZP*S INI I XF-XI ) *DRJ /SI 

Cta F= 1 -CCS ( (XP-XI ) *CR)*CCS( (YP-GAMM AO )*CR ) + S INI I XP-X IJ *DR )*S I NI IYP- 
1GAMM AO)*DR) *CZP )/S I 
V»P = AT AN2I SWF.CWF )*RD 
C 

CALL EULER! >EC, YEC, 2EQ ,XPEQ ,YPEQ ,ZPEQ, 90.+ALPHAO, WP+180. ,90. -GAMMA 

10,0. ,0.,0.,C%,0.,0« ,1,0) 

C 

RE TURN 
END 
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SUBROUTINE RECEQI JC.XEC, YEC.Z EC .XEQ.YEQ.ZEQ ) 

C 

C This SUBROUTINE ROTATES A VECTOR FROM THE MEAN EQUINOX AND ECLIPTIC 

C OF CATE TO THE MEAN EARTH EQUINOX AND ECUATCR OF DATE COORDINATE SYSTEM 

C 

C JC - JULIAN CATE 

C XEC.YEC.ZEC - CCMPCNENTS OF THE VECTOR IN THE MEAN EQUINOX AND ECLIPTIC 

C OF CATE COORDINATE SYSTEM 

C XEC.YEC.ZEC - COMPONENTS CF THE VECTOR IN THE MEAN EARTH EQUINOX AND 

C EQUATOR OF CATE COORDINATE SYSTEM 

C 

REAL JD 

DR=.C17453292519543 
TE=( JD-2415C20.I/36525. 

X I E= 23.452294-0.0 13C12 5*TE-0. OOOOO 164*TE**2-*0.000000503*TE**3 
C = C GS ( X IE *DP ) 

S = S IN { X IE *DP ) 

XEC=XEC 

YEC= YEC*C-ZEC*S 
ZEOYEC*S*Z£C *C 
RETURN 
END 

SUBPCUT INE FXY Z"PQ V ( VX ,VY ,VZ .XI ,W ,0. RPQW, VP. VQ.VW) 

C 

C THIS SUBROUTINE ROTATES A VECTOR FROM THE XYZ TO THE PQM COORDINATE 

C SYSTEM 

C 

C VX.VY.VZ - COMPONENTS CF THE VECTOR IN THE XYZ SYSTEM 

C XI.W.C - INCLINATION, ARGUMENT OF PERIAPSIS, 4N0 LONGITUDE OF 

t THE ASCENDING NODE 

C R PUR - ROTATIONAL MATRIX FROM THE XYZ TO THE PQW SYSTEM 

C VP, VQ.VW - COMPONENTS OF THE VECTOR IN THE PQW SYSTEM 

C 

DIMENSION RPQW ( 3 , 3 ) 

DR=.C174532 '2519943 
C 

CW=CCS IW*CR) 

SW=SI N( W*DR) 

CO = COS ( 0* DR ) 

S0=SIN(0*CRI 
CXI =C0 S t X I *DR I 
SX I=S INIX I*OR ) 

C 

PPCW (1 ,1 >=CW*CO-SW*SO*CXI 
PPCW (1.2) =CN*SC+SW*CC*CXI 
FPCW (1 , 3 ) = SN* SX I 
RPCWI2.1 )=-SW*CC-CW*SO*CXI 
RPC M2, 2) =-SW*SO+CW*CO*CXI 
RPCWI2.3 >=CV*SXI 
RPQW (3,1)=SC*SX I 
PPCW 13 ,2)=-C0*SX I 
RPCWI3 ,3)=CXI 
C 

VP=PPQW<1 ,1 )*VX+RPCW(1 , 2 I *VY+ R PQW ( I ,3)*VZ 
VC=RPQW (2,1 )*VX + RPCW(2,2)*VY«-RPQW(2,3)*VZ 
VW = RPQW(3 ,1 )*V>+RPGW(3,2)*VY+RPQW(3»3) *VZ 
C 

RETURN 
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SUBROUTINE PPCWXYZI VP,VQ,VW, XI » W » 0 * RXY Z » VX ,VY,VZ) 

THIS SLBROUT INE ROTATES A VECTOR FROM THE PCW TO THE XYZ 
COORDINATE SYSTEM 

VP.VOtVW - CCMFCNEMS CF THE VECTOR IN THE FQW SYSTEM 
XI,R,J - INCLINATION, ARGUMENT OF PERIAPSIS, LONGITUDE OF ASCENDING 
NODE 

RXY Z - ROTATICNAL MATRIX FROM THE PQW TO THE XYZ COORDINATE SYSTEM 
VX » VY > VZ - COMPONENTS CF THE VECTOR IN THE XYZ SYSTEM 

DIMENSION RXYZ(3,3 ) 

DR=.C17453292519943 

CR-COS ( W*DR ) 

Sd=SINI W*DR ) 

CO=COS ( 0*DR > 

SO=SIN(C*CR) 

CXI=COS IXI*CR) 

SX I=S I N( XI *CR) 

RXYZ (1 , I )=CV*CO-SR*SO*CXI 

RXYZtl,2)=-SR*CC-CR*SO*CXI 

RXYZI 1,3)=SC*SXI 

RXYZ (2,1 )=CH‘SO + SR*CO*CXI 

RXYZ 12 ,2 )=-SR*SC4CR<C0*CX I 

RXYZI2.3) =-CO*SXI 

RXYZ (3 « 1 ) =SR*SX I 

RXYZI 3,2 ) = CR*S X I 

RXY Z I 3 , 3 ) =CX I 

VX=RXYZ(1 ,1 )*VP+RXYZ( 1 , 2) *VQ+ RXYZI 1 ,3>*VW 
VY = PXYZ 12 , 1)*VP+RXYZ( 2, 2 ) *VQ +R XYZ I 2 , 3 ) *VW 
VZ =R X YZ I 3 , 1 )*VP4RXY2I3 , 2 ) * VO + RX YZ I 3 , 3 >*VW 

RE URN 
END 


SUBPCUT INE EE ARTHI JC,XHE ,YHE ,ZHE, DXHE, DYHE, CZHE) 

THIS SUBROUTINE CCMFUTES THE HELIOCENTRIC POSITION AND VELOCITY OF 
THE EARTH IN MEAN ECUINOX AND ECLIPTIC OF DATE COORDINATE SYSTEM. 
THIS ROUTINE CALLS SUBROUTINES T I NV S AND CCNCAR. 

JD - JUL I AN CATE 

X H E i Y HE ,Z HE - POSITION OF EARTH 
CXH E ,CY HE, CZHE - VELOCITY CF EARTH 

REAL JD 

ANGLE (XI=AMCDI X,3CC . ) +1 8C .-S I GN 1 1 80 . ,X ) 

DR= .074532 92 519 94 3 
RC= 5 7 .2957795130823 
US0N=1.32715445E+11 
AU = 1 49 5 9 8 £45 • 
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C 

0= JD-2415C2C. 

CC=C/1CCCC. 

TE = C/36525 . 

C 

AE = 1 .COOOC023*AL 

Et=C.016751C4-C.00CC418C*TE-0.C00000126*TE**2 
X I E = t . C 

WE=1C 1. 2208 2 2 *C .CCC C47C68*C«-C. 0G00339*CD**2 -*0 .QOC0000 7*CD**3 
CE =C . 0 

XME= ANGLE! 358.475845+C. 9R5600267*D-C.0C 001 12*CD**2 -0.00000 0C7*CD** 
13) 

C 

CALL T INVS(XME*CPtt£»ECE»FE) 

CALL C CNC A R ( At , E E ,X I E ♦ WE ,OE t FE*RO, XEE . YHE , ZE £♦ OXE Et OY HE, OZ EE tl) SUN ) 
C 

RETURN 

ENC 


SLERUIT INE E V ENL S < JC , XH V »YH V ,ZHV , DXHV , DYHV , CZHV ) 

THIS SUBRCL TINE CCMFUTES THE MEAN HELIOCENTRIC POSITION ANO 
VELCCITY OR VENLS IN THE MEAN EARTH EQUINOX ANO ECLIPTIC OF DATE 
CCCFDINATE SYSTEM. T US RCUTINE CALLS SUBROLTINES TINVS AND CONCAR 

C JD - JULIAN LATE 

C XHV.Y HV.ZHV - POSITION OF VENUS 

C DXEV .DYHV »DZHV - VELOCITY OF VENUS 

C 

REAL JD 

ANGL E(X )=AMCD! X.3 6C . ) + 180 .-S I GN I 180. ,X) 

0R=.C17452292519942 
RD=57 .295779513G823 
LSlN=l .32 ?1 5445E+11 
AU=149598845. 

C 

D= JD-2415C2C. 

CD=L/ 1000C . 

TE=D/3£525. 

C 

AV=C .7233216* AU 

E V=C .OC 682 C 69- C . COC 047 7 4* TE 

XI V=2.39363C+0.C01CC5*TE-C.CCCC0097*TE**2 

CV = ? 5. 7 7964 £■*■ C. 8996 50* TE+C . 0C041 1 *TE**2 

HV=120.16C5CC+1.4CEC36*TE-0.000976*TE**2-CV 

XMV= ANGLE (2 12 . 60 3 2 1 9+1 . 602 13C1 540*D+C .( Of 096400*CD**2 ) 

C 

CALL TINVSI >MV*DR»E\»ECV»FV) 

CALL CCNC AR ( A V* EV ,X I V . b V ,0 V ,FV*RD , XHV , YHV , ZFV , DXH V , DYHV , DZ EV ,U SUN ) 
C 

RETURN 

END 
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SUBROUT INE EMAR S( JD ,XHM ,YHM,ZHM, OXHM,DYHM,D 2HM ) 


THIS SUBROUTINE CONFUTES THE MEAN HELICCENTFIC POSITION 
VELOCITY OF MARS IN THE MEAN EARTH EQUINOX AND ECLIPTIC 
COORDINATE SYSTEM. THIS RCliTINE CALLS SUBRCUTINES TINVS 


AND 

OF CATE 
AND CONCAR 


JC - JULIAN DATE 

XHM * YhM , ZHM - POSITION OF MARS 

DXHM.DYHM.DZHM - VELOCITY OF MARS 


PEAL JC 

ANGLE(X)=AHCDIX,360.H-180.-SIGN(180.,X> 

CR=.C17452292519942 

RD= 5 7. 2S577S5 12C 823 

USUN=1.3271E445E*11 

AU=lA9598€45. 


D= JC-2415 02 0 . 
CD=C /I C 000 • 
TE=D/36525. 


AM=1.5236915*AU 

EMM .0533 129C+0. 000 0920 64* TE-O.OC0000077*TE**2 
X IM =1.850 334- C. CCC675* TE+0. 0000 1 2*T E** 2 

GM=4 8.7E6 44i«-C. 77C S « 1+TE-C . CC00015 * TE* *2-0 . C000C 576*TE**3 
WM=334»2182C3+1. 64C759* TE + C .000130*TE **2 -0 . C0000 1 29 *T E**3-0M 
XMM=ANGLE(319 .5 2' 42 5+C. 5240207666*0 +C .00001 3553 *CD**2+0. 00 000002 5* 
1CC**2 I 

CALL TINVSCXMM*CR ,EN ,ECM,FM) 

CALL CCNCARIAM, EM ,X I M , WM , CM , FM*RD ,XHM , YHM , ZPM , UXHM, DYHM, DZ HM ,USUN ) 

RETURN 

END 


SUBROLT INE SUM SX,SV,SZ,SUNB,A,E,XI,W,0,RPCW,U,RS,KK) 

THIS SUBROUTINE COMPUTES THE POSITIONS IN ORBIT WHICH CORRESPOND 
TO VARIOUS LIGHTING ANGLES AND WRITES OUTPUT DATA. THIS ROUTINE 
CALLS SUBROUTINES SUNB AND. CCNCAR, CARSPH, AND TCONIC. 

SX »XY * S Z - COMPONENTS OF THE UNIT VECTOR FROM THE PLANET TOWARD 
THE SUN 

SUNE(l-fc) - LIGHTING ANGLES. (ANGLE BETWEEN THE SUN VECTOR AND THE 
SPACECRAFT RADIUS VECTOR 
A » E - SEMIMAJOR AXIS, ECCENTRICITY 

X I ,H ,0 - INCLINATION, ARGUMENT OF PERIAPSIS, LONGITUDE OF ASCENDING 
NCCE 

RFCW - ROTATIONAL MATRIX FROM THE XYZ TO THE PQW COORDINATE SYSTEM 
U , R S - GRAVITATIONAL CONSTANT AND RACIUS OF THE PLANET 
KK - CONTROL INTEGER. 0 IMPLIES THAT XI,W,0 ARE INPUT. 1 IMPLIES 
THAT RPQW IS INPUT. 


DIMENSION RPQW(3,3) , SUNB ( 6 I 
DR= .C 1 74532S2519943 
ANGRT=SQRT(L/A**3 I 
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MR IT E ( 6 , 1 C5 J 
DC 1C 1=1,6 

IF € SLNB ( I 1 .LT..C11 GC TO ECO 

CALL SUN6ANC ( S X , SY , SZ, SUNB ( 1 1 , XI , W , 0 , RPQW , T A 1 ,T A2 , 1 TYPE1 , I T YPE2 , KK 
1) 

IF I ITYPE1 .EC.C) GC TO 10 

CALL CCNC AR( A,E,XI,k,0,TAl,X,Y,Z,DX,DY,DZ,U) 

CALL CARSPM X, Y.Z.DX.DY ,0 Z ,R 1 ,RA 1, DEC1 , V 1 , F F A l , A Z 1 
CALL T CON IC (U * E , A , T A1 , T1 ) 

Tl=Tl/60. 

H1=R1-R S 

VQHl=Vl*CCS ( FP A 1*GR )/Hl 
GC TC (1,2,3,41 , 1 TY P El 

1 i»R 1 1 E( 6 ,1 Cl ) SUNB (I I ,T 1 , T A1 ,H1 ,DEC1 , RA1 , V0H1 
GO TO 5 

2 WRITE(6,10?) SLNB(I ) ,T1,TA1,H1,DEC1,RA1,V0H1 

GO TC 5 

3 InRI TE 16 ,103) SUNB(I),T1,TA1,H1,DEC1,RA1,V0H1 


GO TO 5 

4 URITE(6,104) SUNB ( I ) , T 1, T A 1 ,H 1 ,DEC 1 , R A 1 , VOH 1 

5 CALL CCNC«R(A,E,XI,lrf,0,TA2,X,Y,Z,DX ,DY ,DZ,U 1 

CALL CARSPH(X,Y ,Z ,CX,DY,DZ,R2,RA2,DEC2 ,V2,FFA2,AZ) 
CALL TCCNIC1U ,E ,A,T A2,T2) 

T 2 = T 2/6 0 . 

H2=R2-RS 

V0F2=V2*CCS(FPA2*DR)/H2 
GO TO (6,7,8,91, ITYFE2 

6 WRITE(6,101) SUNBd 1 ,T2,TA2,H2,0EC2 ,RA2,V0H2 
GO TC 1C 

7 WRITE(6,102) $LN3(I),T2,TA2,H2,DEC2,RA2,VOH2 
GO TO 10 

8 WR IT E ( 6 ,1 C3 1 SU NB ( I) , T2 ,T A2 ,H2 , OEC 2, R A 2, V CH2 
GO TO 10 

9 WR IT E ( 6 , 1 C4 1 SUNB (I 1 ,T2 ,TA2 ,H 2 , DEC 2 , RA2 , V0H2 


10 CONTINUE 


101 

FOFMAT ( 11 H 

SUN 

ANGLEF5.1.21H 

DEG 

AND 

INC , 

s/c 

ASC,5Fll.2,F14.8) 

102 

FORMAT ( 11 H 

SUN 

ANGLEF5. 1, 21H 

DEG 

AND 

INC, 

S/C 

CSC,5F11*2.F14.8) 

103 

FORM AT ( 1 1 H 

SUN 

ANG LEF 5 • 1 » 2 1H 

□EG 

AND 

DEC, 

S/C 

A SC t5 FI 1*2 »F14*8) 

104 

FORM AT ( 11 H 

SUN 

ANGLEF5.1,21H 

□ EG 

AND 

DEC, 

S/C 

DSC,5F11.2tF14. 8) 

105 

FORMAT ( LHO , 

41X, 

6 2 FT I ME T 

• A* 


ALT 

DEC 


1 V/H/l 


RA 


600 RETURN 
END 
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SUBROUTINE SUNe ANC ( S X , S Y , SZ , P S I , X I ,W,0,RPQW ,T A 1 ,T A2 , IT YPE1 , ITYPE2, 
IKK) 

C 

C THIS SUBROUTINE CCMFUTES TEE POSITION IN OR El T WHICH CORRESPONDS 

C TO A GIVEN LIGHTING ANGLE. THIS ROUTINE CALLS SUBROUTINES RXYZPQW 

C ANC CACRAT 

C 

( SX » S Y 9 SZ - CCMFCNEMS OF THE UNIT VECTOR FR CM THE PLANET TOWARDS 

C THE SUN 

C PSI - LIGHTING ANGLE. (ANGLE BETWEEN THE SUN VECTOR AND THE 

C SPACECRAFT RACIUS VECTOR) 

C XltWtO - INCLINATION, ARGUMENT OF PERIAPSIS, LONGITUDE OF 

C ASCENDING NODE 

C RPCW - ROTATIONAL MATRIX FRCM THE XYZ TO THE PQW COORDINATE SYSTEM 

C T A 1 * T A2 - FIRST AND SECOND TRUE ANOMALIES C CRR ESPOND I NG TO A 

C LIGHTING ANGLE 

C ITYPEI , ITYPE2 - CCNTROL INTEGER TO INDICATE THE CONDITIONS AT T A 1 

C ANO TA2. 0 IMPLIES NO SOLUTION. I IMPLIES PSI INCREASING 

C SPACECRAFT ASCENDING. 2 IMPLIES PSI INCREASING, S/C 

C CESCENDING. 3 IMPLIES PSI DECREASING, S/C ASCENDING. 

C A IMPLIES PSI CECREASING, S/C DESCENDING. 

C KK - CO NT PC L INTEGER. 0 IMPLIES THAT XI,W,0 ARE INPUT AND RPQW IS 

C OUTPUT. 1 IMPLIES THAT RPQW IS INPUT. 

C 

DINENSICN F F C W ( 3 ,3) 

CR = .0 17 A5 32S2 5 1 ^5A 3 

I T V FE 1 =0 

ITYPE2=0 

IF (KK «EQ • 1 ) GO TO 1 

CALL RXYZPCMSX ,SY,SZ,XI ,W ,0 ,RPQW ♦ SP, S Q, S W) 

GO TO 2 

1 SP=RPQM1 , 1J*SX+RPCW(1 ,2)*SY+RPQW(I ,3) *SZ 
SQ = R PQW (2,1 l*SX+RPCW (2 ,2) *S Y + RPQW ( 2 , 3 ) * SZ 
SW = PPQW(3 ,1 )*SX+RPCW(3, 2)*SY+RPQW(3,3)*SZ 

2 A=SQ**4/SP**2+2 .*SQ**2+SP**2 

B = -2 .*CCS(PSI *DR)*(SQ**3/SF**2+SQ) 

C=SW**2*SQ**2/SP**2*C0S( PSI*DR)**2-SIN(PS I*CR )**2 
CALL QADRAT(A,B,C,SF1,SF2,NUM) 

IF I MM.EQ.C )GC TO 800 

IF( AESISFi ) .GT .1 ..CP. ABS(SF2> .GT . 1)G0 TO EC C 
TA1=ATAN2(SF1 , (COS( PSI*CR )-SC*SFl )/SP)/DR 
T A2 = AT AN2 ( SF2 , (CC$( PSI*DR )-SQ*SF2)/SP)/DR 
SOPS I=( SP*SF1-SC*CCS(T A1*CR) )/ S I N ( PS I*DR ) 

A SC E S =RPQW (2,3)*CCS(TA1*DF)-RPQW(1 ,3)*SF1 
IF l SDPSI.GT.O. . ANC. ASDES.GT .0. )I TYPE1 = 1 
IF ( SCPSI.GT.O. . ANC. *SDES.LT .0. ) ITYPE1=2 
IF(S£PSI.LT.O..ANC.ASDES.GT.O. )ITYPE1=3 
IF(SDPSI.LT.O.. ANC. ASDES.LT .0. )ITYPE1=A 
SOPS I= t (SP*SF2-SQ : *CCS(TA2*CR) )/SIN(PSI*DR ) 

ASDE S = FPQ W ( 2 , 3 ) *CCS ( TA 2*D R )- RPQ W (1 ,3) *SF2 
IF I SDPSI .GT.O. . ANC. ASDES.GT .0. ) ITYPE2=1 
IF (SCPSI.GT.O. .AND. ASDES.LT. C. )I TYPE 2= 2 
IF I 5 CP SI .LT.O. .ANC. ASDES.GT .0. )ITYPE2=3 
IF(SCPSI .LT.O.. AND. ASDES.LT. 0. )ITYPE2=4 
800 CONTINUE 
RETURN 
END 
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SUEROUTINE OCC ULT ( A , E , XI , W , 0, U, RS , EX , EY ,EZ ,TOCC ,T 1 , ALT1 ,F 1 , DEC1 , R 
1A1 ,T2,ALT2,F2,DEC2,RA2,KK) 

C 

C THIS SLBROLT INE COMPUTES THE ENTRANCE AND EXIT TRUE ANOMALIES OF 

C OCCILTATICN. THIS ROUTINE CALLS SUBROUTINES RXYZPCW, QARTIC, 

( CROSS, DOT, TCCNIC, RPQWXYZ, AND LATLNG * 

C 

C A , E , X I - SEMIMAJOR AXIS, ECCENTRICITY, INCLINATION 

C W,0 - ARGUMENT OF PERIAPSIS, LONGITUDE OF ASCENDING NODE 

C U » RS - GRAVITATIONAL CONSTANT AND RADIUS OF THE PLANET 

C EX»EY»EZ - COMPONENTS OF UNIT VECTOR TOWARD THE BODY OCCULTED 

C TOCC - LENGTH OF TIME IN SHADOW 

C T1»ALT1»F1,CEC1»RA1 - CONDITIONS AT ENTRY INTO THE SHADOW, TIME 

C FRCM PERIAPSIS, ALTITUDE, TRUE ANOMALY , DECLINATION, AND 

C RIGHT ASCENSION 

C T2»ALT2»F2»CEC2,RA2 - CONDITIONS AT EXIT FRCM THE SHADOW 

C KK - CONTROL INTEGER. 0 IMPLIES NO OCCULTATION, 1 IMPLIES OCCULT. 

C 

DIMENSION PPUW! 3,3) ,CF!4),RXYZ(3,3> 

ANGLE! X )=AMCD( X »36C . ) + 18C .-S I GN ( 180. ,X) 

DR = .C17452292519 < sA2 
RD=57.2S577«5 13C823 
F1 = 2CCC • 

F2 = 200C . 

K K = 0 

P = A* ( 1- E*E ) 

CALL RXYZPJMEX,EY,EZ,XI , W , 0 , RPQW , BET A ,XX I ,ZBODY) 

Cl=(RS/P)**4*E**4-2.*( RS/ P ) *#2* ( XXI **2-BETA**2 >* E**2 M BETA **2 +XX I * 
1 * 2 ) * * 2 

C2=4. *(RS/P)**4*E**2-4.*( RS/ P ) **2* ! XXI **2-E ETA** 2 )*E 
C3=6.*(RS/P)**4*E**2-2.*(RS/P)**2*( XXI **2-BETA**2)-2.*(RS/P)**2*( 1 
l.-XX l**2)*E**2*2.*( XXI**2-BET4**2)*( 1 .-XX I * * 2 ) -4. *BE TA**2* XX 1**2 
C4=4.*(RS/P )**4*E-4 ,*(RS/P)**2*( l.-XXI**2)*E 
C5=(RS/P) **4-2.*(RS/P)**2*( 1 .-XX 1**2 > + ( 1 . -X XI **2 ) **2 
CALL CART IC (Cl,C2,C3,C4,C5,CF(l ) ,CF (2 ) ,CF( 3 ), CF( 4 ), JJ ) 

IF(JJ.LG.O) GO TO EOC 
C 

CO 3 1=1, JJ 

IF ( ABS( CF ( I ) I.LT.1.C001. ANC.ABSICF! I)) .GT.l .) CF( I )=0.999999 
IF(AI3S(CF( I H.GT.l. ) GC TO 3 
SF = SCRT ( l.-CF! I )**2 ) 

R = P/(1 . +E *CF ( I ) > 

CALL CROSS ( R*CF 1 1 ) ,R*SF»0 «,BET A, XXI , ZBODY.P >, PY, PZ .PRODUCT ) 

CALL DOTCCF! I ) ,SF ,C . .BETA ,XX I , Z BODY , ANG ) 

I F ( A ES ( PRO DUCT- RS ) .IT. .Cl. AND. ANG.GT.90. ) GC TO 1 
SF=-SF 

CALL CROSS (P*CF C I ) , P *S F ,0 ., BET A , XX I , ZB OD Y , P >, PY, PZ, PRODUCT ) 

CALL DCT(CF(I),SF,0., BETA, XXI ,ZBODY , ANG) 

IF(ABS(PRODLCT-RS).LT. .01. AND. ANG. GT .9C.) GO TO 1 

GO TO 3 

1 IFIFl.LT. 10 CO.) GO TO 2 

Fl= ANGLE! ATAN2 (SF,CF(I ) )*RD) 

GO TO 3 

2 F2=ANGLE< ATAN2(SF,CF( I ) )*RD ) 

GO TC 4 

3 CONTINUE 
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4 IF<F2.GT.1CC0.)GC TC 800 
CF 1= COS ( F 1*CR ) 

SF 1= S IN < F 1 ♦DR ) 

DS CF =2 • *R S*R S * ( 1 . + E *C F 1 >* l-E*SF 1) + 2 . *P* P* C B ET A*CF 1 + XX 1 * SF1 ) * < -BE T A 
1*SF1+XXI*CF1 ) 

IFiDSDF.GT.C. )G0 TO 5 
C 

T EMP = F 1 
F 1 =F2 
F2-T EMP 
C 

5 CALL TCCNICCU.E ,A,F1,T1I 
CALL ICON IC(U»E«A»F2vT2) 

T CCC= T 2-T 1 

IFiTOCC .LT.C . )TCCC=6.28 3185 3*SQRT( A**3/UH-TCCC 
T 1 = T 1/ 60 . 

T2=T2/60. 

T0CC=T0CC/60. 

AL 1 1 = P/ ( 1,+E*CCS(F1*DR) )- PS 
ALT2=P/ ( 1 • + E* COS ( F 2* DR ) )-RS 

CALL RPGWXYZ(CCS( FI ♦DR ) ,S IN(F1*DR) , 0 . , X I * W # C, RXYZ , RX , R Y ,RZ ) 

CALL LATLNG(RX,RY,RZt DEC1,RA1 ) 

CALL RPCWXYZl CCS(F2*CR ) ,S IN ( F 2*DR )> 0 . , XI ,W, C, RXYZ,RX ,RY ,RZ ) 

CALL LATLNG (RX, RY,RZ,DEC2 ,RA2 ) 

KK— 1 

GO TC SCO 
C 

ECO CONTINUE 


T 0 C C =0 . 
T1 = C • 

AL T 1 =C • 

F 1 =C • 
DEC1 =C • 

R A 1 = 0 • 

T2 -C • 

A L T2 = C • 
F2 =0 • 
DEC2 = 0 • 

R A2-C • 

900 CONTINUE 
RETURN 
ENC 
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SLERULTINE CONC AR ( A , E , X I , W ,0,F , X, Y. Z , DX , DY » EZ ,U ) 

C 

C THIS SLBROOT INE CONVERTS CCNIC ELEMENTS TO CARTESIAN COORDINATES 

( 

C A # E t X I - SEMMAJOR AXIS, ECCENTRICITY, INCLINATION 

C UfO'F - ARGUMENT CP PERIAFSIS, LONGITUDE OF ASCENDING NODE, TRUE 

C ANLMLY 

C X » Y , Z - CCNPCNENTS CF POSITION VECTOR 

C CX fCY tDZ *■ CC RFC NE M S OF VELOCITY VECTOR 

C U - GRAVITATIONAL CONSTANT 

C 

CATA CR/ *017453292519943/ 

FR=DR*F 
kFR=UR*(W4F ) 

CR=DR*G 

XI F=DP*XI 

C E N- 1 • *E*CC 5 ( FF ) 

R = A*< 1 .-E*E )/DEN 
V=SCRT(U*<2 ./R-l./ A )) 

GAM=ATAN( E* S I N( F P ) /DEN ) 
kFGP^kF R-G AN 
CWF=CCS< UFRJ 
ShF = SIN(WFR ) 

SC= S I N ( CR) 

CO=CC S ( OR ) 

S I = S IN ( X IR ) 

C I =COS ( XI R ) 

S1nFC = S I N ( WFGR ) 

CKFG=COS< UPGR) 

X=P* <CWF*COSkF*SO*C I) 

Y=R 5 MCWF*SC4 SWF *CC *C I ) 

Z=R*SkF*S I 

DX=V*(-$WFG*CC-CWFG*SO*CI » 

OY=V*(-$wFG*S0+CWFG*CO*Cl ) 

D Z= V*CWFG*S I 

RETLRN 

END 
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SUBPCLT INE CARSPHIX ,Y ,Z ,DX,DY,DZ,R, R AS , DEC , V , G AM , S IG ) 

C 

C THIS SUBROUTINE CONVERTS CARTESIAN COORDINATES TO SPHERICAL 

C CCCPDINATE. 

C 

C X » Y » Z - COMPONENTS CF PCSITICN VECTOR 

C OX.CY.DZ - COMPONENTS OF VELOCITY VECTOR 

C R?RIS *CEC - MAGNITUDE CF POSITION VECTOR, RIGHT ASCENSION, 

C CECLINATION 

C V.CAM.SIG - MAGNITLCE OF VELOCITY, FLIGHT-PATH ANGLE, AZIMUTH 

C 

DATA RD/57. 2 S57 7 S5 1 208232 1 / 

ARSM T)=ATANU/SCRT (1 .-T*T» ) 

R=SQRT(X*X+Y*Y+Z*Z ) 

R A S= AT A N2 I Y • X ) 

DEC= ARSNl Z/R ) 

V= SCRT (DX*CX+DY*CY+£Z*DZ ) 

CT =CC S ( R A S I 
ST =S INI HAS ) 

CP=C OS (DEC ) 

SP=S INI DEC ) 

DXP=CX*CT*CP+DY*ST *CP+DZ*SP 

DYP=-DX*ST+CY*CT 

CZP = -CX*SP*CT— D Y*SP<ST+CZ*CP 

GAM=RC*ARSN(CXP/V) 

SIG=PD*AT AN 2 I CYP , DZP) 

R AS=RC*RAS 
C£C=RC*DEC 
RETUPN 
END 
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SU8PCUT INE TCONIC(U.EC,A,TA22,T) 

C 

C THIS SUBROUTINE COMPUTES TEE TIME FROM PERIAPSIS PASSAGE FCR A 

C GIVEN TRUE ANCMALY 

C 

C U - GRAVITATIONAL constant 

C EC.A.TA22 - ECCENTRICITY, SEMIMAJOR AXIS, TRUE ANOMALY 

C T - TIME FPCM PEPIAFSIS PASSAGE 

C 

T AN GF ( X )= S l N (XI/CCS IX) 

TA2=T422*.174S2 2<;2 5E-1 
SLR= A*( I * -EC *EC ) 

AE= ABS (A) 

FAC= AB^SORT IAB/U) 

EC A= ( 1 .-EC )/ll . +E C ) 

ABE = SORT < AES ( ECA) ) 

THt = T ANGF I.5*TA2) 

IFIAEE-.CCCC5 111 ,11 ,12 

12 CONTINUE 

EC A=2 . *AT AN (AEE*TF£) 

I F ( A) 14,11 ,13 

13 T=FAC*< ECA-EC*S IN (ECA)) 

GO TC 16 

14 ANG = ABE *T FE 
ANG=1.+2.*ANG/ ( l.-ANG) 

T = F A C *( EC *TANGF ( EC A )-ALOG( ANG) ) 

GO TO 16 

11 F AC= SQR T <SLR**2/U)*2./< ( l.+EC )**2) 

EC 1 = EC A *THE * *2 

T=F AC*( THE + TFE**3M (l.-2.*ECA)/3.-(2.-3.*ECA)*ECl/5.+(3.-4.*ECA>*E 
1C 1**2/ 7.- (4 .-5.* EC A )*ECl**3/<5. ) ) 

16 CONTINUE 
RETURN 
END 

SUBROUTINE T IN\ SIM ,E , EC • F ) 

C 

C THIS SUEROLTINE CONVERTS MEAN ANCMALY TO ECCENTRIC AND TRUE 

C ANCM4LY 

C 

C K , E - MEAN ANCMALY AND ECCENTRICITY 

C EC , F - ECCtNTRIC ANO TRUE ANOMALY 

C 

REAL M , MU 
EC =M 

10 MO-=EC-E*S IN(EC) 

DM=M-MG 

Dfc-CM/ ( 1 • — E*CC $ I EC ) ) 

EC = EC +DE 

IF( ABSIOE ).GT. .0001 )G0 TC 10 
H EC = EC/2. 

HF = A T AN ( SCR T CU. + E)/C1.-E) ) *S I N( HEC ) /CUS I HE C ) ) 

IFIHF.LT.C.) HF =3. 1415926+HF 

F=2.*FF 

RETURN 

END 
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SUBRCUT INE CACR «TC« ,B,C ,X I ,X2,KK) 

C 

C THIS SUBROUTINE SOLVES THE EQUATION AX**2 + EX +C=0 FOR THE REAL 

C ROOTS 

C 

C A,B,C - COEFFICIENTS OF TFE DIFFERENT POWERS OF X 

C XI, X2 - REAL ROOT OF THE EQUATION 

C KK - NUMBER OF REAL ROOTS 

C 

KK= C 

CIS = B*B-4.*A*C 
IFICIS.LT .0 .) GO TO 800 
XI a I -B ♦SORT (DISII/2./A 
X2=( -B-SQRTCD IS ) )/2 ./A 
KK = 2 

£00 RETURN 
ENC 

SUBROUTINE CUEICIA.B.C.D, XI, X2 , X3 » K K ) 

C 

C THIS SUBROUTINE SOLVES THE EQUATION AX**3 *EX**2 +CX +D = C FOR 

C THE REAL ROCTS 

C 

C A,E,C,D - COEFFICIENT OF THE DIFFERENT POWERS OF X 

C X1.X2.X3 - PEAL PCCTS OF THE EQUATION 

C KK - NUMBER CF REAL ROOTS 

C 

CBRT ( X ) = S IGMAB SIX)**. 2332-323 3, X) 

KK=0 

PI=3. 1A15927 

IFIA.LT .. 1E-30.ANC. A. GT.-. IE-301 GO TO A 
P=B/A 
C=C / A 
R = D/A 

SA=(3.*Q— P**2)/2. 

SB=(2.*P**3-9.*P*Q+27.*R)/27. 

DEL = ( A .*Q **3-Q**2*F**2-18.*Q*P’*R*2 7.*R**2+A .*P**3 *R ) /108. 
IF(CEL.LT..lE-30.AND.DEL.GT.-.lE-30) GO TO 2 
IF ( DEL 11,3,2 

1 KK= 3 

CPHI=-SB/2./SQRT ( S A**3 / 1 -27 . ) I 
IFIABSICPHI J.GT . 1 . ) C-0 TO 1C 
SPHI=SQRT ( 1 .-CFHI**2 I 
PHI = AT AN2(SPHI,CPHI ) 

GC TC 11 

10 SPHI=SCRT((27.*DEL)/SA**3) 

C SINCE FOR SPALL ANGLES SPHI=PHI 

BETA=SPHI 

IF I-SB .GT .0 • ) PH I = BE T A 

IF l - SB. LT.O.) PH I =3.14159 26 5 3 589793-BET A 

11 E0 = 2 •* SQR Tl-SA/3.) 

XI = EO*COS I PH 1/ 3 . )-F/3. 

X2=EC*CCS(FHl/3.+2.*PI/3.)-P/3. 

X3=EO*COS(PHI/3 .♦4.<PI/3.)-P/3. 

GO" TO 7 

2 KK= 1 

X 1 =C ERT (- SB /2 . + SCRT (DEL) ) +CBRT (-SB/ 2.-SQRT CCEL) I-P/3. 

GO TO 7 

3 KK=3 

X 1 = 2. *CBRT( -SB/2. l-P/3. 
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X2=C6RT< SB/2 . )-P/3 • 

X3 = X 2 
GC 1C 7 

4 CONTINUE 

DI S=C**2-4.*B*D 
IF(0ISi7,5,5 

5 X 1 -4-C +SQRT ( D I S ))/2*/8 
>2=(-C-SGRT(DIS))/2./B 
KK=2 

7 CONTINUE 
RETURN 
ENC 


SUBROUTINE Q ART 1C( A ,B,C ,D ,E , XI , X2 , X3 , X4 ,KK) 

C 

THIS SUBROUTINE SOLVES TEE EQUATION AX**4 *EX**3 *CX**2 +DX *E = 0 
FOR THE REAL ROOTS. THIS ROUTINE CALLS SUBROUTINES Q ADR AT AND CONIC 

A,fitC,G,E - COEFFICIENTS OF THE DIFFERENT POWERS OF X 
X1»X2«X3«X4 - REAL ROOTS OF THE EQUATION 
KK - NUMBER OF REAL ROOTS 

KK = C 
B P-E/ A 
CP=C/A 
CP = C / A 
E F = E / A 
C 

P/4 • 

H2 = H **2 
F3 = F 2 * H 
H4=F3*H 

F=6.*H2*3.*EP*F+CP 
Q=4.*F343 # *BP*H2*2 .4CP*H+CP 
R=H4 +BP*H3+CP*F2*DP*H + EP 
C 

CALL CUB IC< 1 .»2 .*F, R *P-4 • *R ,-Q*Q ,T1 ,T2»T3, NFCOT ) 

C 

GO TC (1,2,3) , N ROOT 

1 RP =T 1 
GO TO 4 

2 RP = A MA Xl ( T1 * T2 ) 

GO TO 4 

3 RP= AMAXK T1 ,T2 , T3) 

C 

4 CONTINUE 
SCFP^SGKT (RP) 

XI MP+RP-G/SQRP )/2 . 

BETA=(P+RP+C/SQRP)/2. 

C 

CALL UACR AT ( 1 . , SCR P , X I , Y 1 , V 2 , 1 ROOT ) 

CALL QACRATCl. ,-SGP P ,BE TA , Y3 , Y4 , JROOT ) 

IF ( IRCCT+ JRCC T • EC .0 ) GO TO 800 
I F 1 IROGT + JRCCT .EC. 4 ) GO TO 6 
I F ( IROOT.EQ.O ) GO TC 5 
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X I=Y 1+H 
X2=Y2+H 
KK=2 

GO 10 800 

5 X 1= V3+H 
X2 =>4+H 
KK=2 

GO TO 800 

6 X 1 =Y1 
X2 =Y2*H 
X3=Y3+H 
X4=Y4*H 
KK=4 

€00 CONTINUE 
RETURN 
ENC 


r 

SLBRCUT INE 

CRGSSCXl » Y 

ItZ 

1 *X2 

t Y2f Z2 » 

PX* PY,PZ , PRODUCT) 

C 

r 

THIS SUBROUTINE COMPUTES 

THE 

VECTOR 

OR CRCS5 PRODUCT 

V 

c 

X1,Y1,Z1 ~ 

CCMPCNEMS 

OF 

THE 

FIRST 

VECTOR 

c 

X 2 * Y2 • Z 2 - 

CC ► FCfvEN 1 S 

OF 

THE 

SECOND 

VECTOR 

c 

PXtPYtPZ - 

CCPPCNEMS 

CF 

THE 

VECTOR 

PRODUCT 


C PPCCUCT - MAGNITUDE OF THE VECTOR PROOUCT 


C 

PX =Y 1*Z 2-Z 1*Y2 
PY=Z1*X2-X1*Z2 
PZ=X 1*Y2-Y1*X2 

PPCDUCT=SQRT( PX*FX+FY*PYtFZ*PZ) 

RETURN 

ENC 


SUBROJT INE COT (XI ,Y I ,Z 1 ,X2 ,Y2t Z 2, ANGLE ) 

C 

C THIS SUBROUTINE COMPUTES THE ANGLE BETWEEN TWO VECTORS 

C 

C XltYltZl - COMPONENTS OF THE VECTOR R1 

C X2.Y2tZ2 - COMPONENTS OF THE VECTOR R2 

C ANGLE - ANGLE BETWEEN VECTORS RI AND R2 

C 

R0=ET.2S5TT<513Ce23 

R1=SQRT(X1*X1+YI*Y1+Z1*Z1» 

R2=iCPT (X2*>2+Y2*Y2*Z2*Z2 ) 

C 

ANGL E = ACOS ( I X 1*X2-»Y 1 =»Y2*Z1 * Z2I /R1/R2I*RD 

RETURN 

ENC 
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SLEROUT INE CRB I T ( SL AT , SLGNG , VI NF ,HA ,HP , BET A , A , E , XI ,W,0,PLAT ,PLONG, 
1U.RS ) 

THIS SUBROUTINE CONFUTES TEE KEPLER IAN ORBITAL ELEMENTS FOR A 
HYPERBOLIC PERIAPSIS TO AN ELLIPTICAL PERIAPSIS TRANSFER 

SLAT, SLUNG - LATITUDE AND LONGITUDE OF THE S-VECTCR 
VINF - HYPEPEOLIC EXCESS VELOCITY 

FA, HP - AFCAPSIS ANC PERIAPSIS ALTITUDES OF THE ELLIPSE 

BETA - ANGULAR MEASUREMENT OF THE ORIENTATION OF THE ORBITAL PLANE 

A , E ,X I - SEMINA JO R AX I S , E CC ENTR IC I T Y , AND INCLINATION 

H,G - ARGUMENT OF PERIAPSIS AND LONGITUDE OF THE ASCENDING NODE 

PLAT »PLCNG - LATITUCE ANC LONGITUDE OF PERIAPSIS 

U , RS - GRAVITATIONAL CONSTANT AND RADIUS OF THE PLANET 

ARCOS(X) = ACCS( X J 
ARSINIX ) = AS IN ( X ) 

ANGLE! X )=AMCD (X ,360. H- 18C ,-S I GN( 1 80 . , X ) 

RC=57 .255775513C823 
DR =.C 17 45 2 2 525 15 542 
X I=BETA 

I F (X I .GT. 85 .5 555 .AND. XI .LT. 50. 0001) XI =85.5 599 
IF (XI. GT. 265.5559.AND.XI.LT. 2 70. CCC1) XI =2 6 5. 9999 
A= ( 2 .*PS+FA-*HP ) /2. 

E= ( R S + HA ) /A-l . 

APSHT = ARCGS IU/ (U+(RS+HP )*VINF**2 ) )*RD 
B=ARSIN(SIN(SLA1*DP>/AES(SIN(X I*DR) ) )*RD 
DEL = ARS IN< TANI SLAT*CR )/ABS ( TAN (XI* DR) ) )*RD 
IF (XI. LT. 180.) GC TO 2 
W=AN GLE ( 1 6C.-B- APSHT) 

IF (XI .LT. 270. ) GC TC 1 
0= ANGLE! 180.+SLCNG+CEL ) 

GO TO 4 

1 0= AN GLE (18C.+S LONG-CEL) 

GO 10 4 

2 N* ANGLE ( B-APSHT ) 

IF ( XI . L T. 90 . ) GO TO 3 
0= AN GLE ( SLCNG + DEL ) 

GO TO 4 

3 C= ANGLE (SLCNG-C EL) 

4 CONT INCE 

IF(X 1 .GE.18C. ) X 1= 2 6C. - XI 

PL AT=ARSIN ( S I N ( W*CR )*S INI X I *DR ) ) *RD 
SSDEL = T AN ( PL AT* CR ) / T AN ( X I*DR ) 

CSOFL = SSOEL/TANU*CF)/CCS(XI*DR) 

SDEL = AT AN2( SSCEL,CSCEL )*RD 
PLCNG= ANGLE (O+SDEL) 

RETURN 

END 
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SUBROUTINE PR EC £S < JC1 ,XE1 t YE1, ZE1 , J02 , XE2 , YE2 , ZE2 ) 

C 

C THIS SUBROUTINE TRANSFORMS MEAN EARTH EQUINCX ANO EQUATOR COORDINATES 

C FROM CNE EPCCH TO ANOTHER EPCCH 

C 

C JC1.JD2 - JULIAN CATES CF INITIAL ANO FINAL EPOCH 

C X E 1 tVEl * Z E 1 - CCMPCNENTS CF VECTOR IN JDl COORDINATE SYSTEM 

C XE2 ? YE2 tl E2 - CCMPCNENTS CF VECTOR IN JD2 COORDINATE SYSTEM 

C 

REAL JC 1 » JD2 

T^ABSl ( JC2-JD1 )/3E524. 219879) 

TO=< JD2-241502C.)/3*524.219879 

ZETAC=< 0. 64CC6944+C .38 7?7778E-3*T0 ) *T+C.83E E8889E-4*T**2+0 . 5E-5*T* 

1*3 

CZETAO-ZETAC+C .21972222E-3*T**2 

THETAC=(0 .556856 11-C. 2 3*9 44 4E- 3*TO ) *T-C. 118 23333E-3*T**2-0 .1166666 
17E-4 *T **3 
C 

IF ( JC2-JD1.GT.0.) GC TO 1 
TE MF= ZE T AO 
ZETAO=-CZET AC 
CZET AC = -T EM P 
THETAO=-THETAO 
C 

1 CALL EULER <> El • V E 1 * Z El » X E2 , YE2 , ZE 2 , 9C . -ZET A ( # - ( 90 . +CZ E T AO ) ? THET AO » 
lDPHtDP SI .DPSI* HXP ,bYP »WZPt 1.01 
C 

RETURN 

END 


SUBROUTINE L A TL NG ( X , Y * Z , X L AT f XLCNG ) 

C 

C THIS SUBRCL TINE CCMFUTES THE LATITUDE AND LONGITUDE OF A GIVEN 

( POSITION VECTOR 

C 

C X,V,Z - CLMFCNENTS CF THE POSITION VECTOR 

C XL AT tXLON G - LATITUDE AND LONGITUDE 

C 

APCCS 1 X ) = ACC S ( X ) 

AR S IN ( > ) = A S IN IX ) 

RD=57. 2^57795 
R=SCRT(X**2+Y**2+Z**2) 

XLCNG = ATAN2 (Y,X )*PD 
XLA7=ARSIN (2/R)*RC 
RETURN 
ENC 
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SOEsRQUT INE EULER (X , Y, Z ,XP ,YP ,ZP , PH I , PS I , THE 1 A ,DPHI ,DPSI »DT H ET A »W XP 
1 ,V*YF ,WZP, J, * ) 

C 

C THIS SUBRCL TINE PEPFCRMS AN EULER ROTATION AND/OR COMPUTES THE 

C INSTANTANEOUS ANGULAR VELOCITY VECTOR 

C 

C X t Y * Z - COMPONENT S CF THE VECTOR IN THE UNP P I MED COORDINATE SYSTEM 

C XP , YP » ZP - COMPONENTS OF THE VECTOR IN THE PRIMED COORDINATE SYSTEM 

C PHI, PSI .THETA - THE THREE EULER ANGLES 

C CPH I.CPSI.CTHETA - TIME RATE OF CHANGE OF THE EULER ANGLES 

C RXP.HYP.UZP - THE PRIMED COMPONENTS OF THE INSTANTANEOUS ANGULAR 

C VELOCITY VECTOR 

C J - CCNTRCL INTEGEP. -1 IMPLIES XP.YP.ZP INPUT, X.Y.Z OUTPUT. 

C 0 IMPLIES NC RCTATICN. I IMPLIES X.Y.P INPUT, XP.YP, ZP CUT PUT . 

C K - CCNTRCL INTEGER. -1 IMPLIES WXP.WYP.WZP INPUT, DPHI.DPSI, 

C DTHETA OUTPUT. 0 IMPLIES NO ANGULAR VELCCITY COMPUTATIONS. 

C 1 IMPLIES CPHI, CPS I, DTHETA INPUT, NXP,V*YP,WZP OUTPUT. 

C 

XPFI=PH1*. 0174532925 
XPSI=PS I*. C 174522 925 
XTH=THETA *.01745329 25 
I F ( J ) 1 C ,12,11 

10 X= (COST XPSI )*CCS(XPHI )-CCS( XTH)*S INI XPHI * * S IN C XPSI) I *XP ♦ (- S I N( XPS I 
1 ) *C C S I XPHl I-CCS (XTH)*SIN( XPHI )*COS( XPS I I) * Y F* ( SI N ( XTH) *S IN C XPH I ) ) * 
2ZP 

Y= (CCS( XPS I )*S IN (XPHI I + COS( XTH I *COS ( XPH I )*SIN ( XPSI ) )*XP+(-SIN( XPSI 
1 >*SIN( XPHI I4CCS ( XTH I *CCS (XPHI )*COS(XPS 1 1 ) *Y F+ ( -S IN ( X TH )*COS ( XPHI ) ) 
2*ZP 

Z=< SIN( XTH) *SIN(XPS I) )*XP+(SIN(XTH)*COS( XPS 1) ) *YP+(COS(XTH ))*ZP 
GO TO 12 

11 XP= (COS (XPSI J*COSI XPHI )-CCS( XTH ) *SI N ( X PH I > * S I N ( XPS I > I *X + (C C S < XPS I I 
1*S IN ( XPHI ) + COS( XTH) *COS( XPHI >*SIN(XPSI ) ) *Y + < S I N< X TH ) *S IN (X PS I)) *Z 

YP= l-S IN ( XPSI >*C0S( XPHI )-COS(XTH)*SIN( XPHI ) *C OS ( XPS I ) ) *X* ( - $ I N< X PS 
1 11*5 IN ( XPHI )+CO SIXTH )*CCS( XPHI )*CCS<XPSI ) ) * U ( COS ( XPS I >*S I M XTH ) )* 

21 

ZP=(SIN(XTH)*SJN(XPHI))*X+(-SIN(XTH)*COS(XPHI) ) * Y+CO S( XTH) *Z 

12 IFIK)13,15,14 

13 CPH l=(WXP*SIN(XPSII*WYP*CCS(XPSm/SIN( XTH) 

DPSI=kZP-(CCS( XTH )*(WXP*SIMXPSI ) +WYP*COS(XFS I) ) )/SIN(XTH) 
DTHETA=fcXP*COS (XPSI )-WYP*SIN( XPSI ) 

GO TC 15 

14 WXP = CPFI*SIM XTH MS INI XPS IT* DTHETA* COS (XPS I ) 

NYP=CPHI*SIM XTH.MCC S( XPSI ) -DTHET A*S IN ( X PS I ) 

WZP=DPSI+DPHI*COS (XTH) 

15 RETURN 
ENC 
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SIBROUT INE VECTCR1 JD,OEC S , R AS ,OECE ,R AEt DECC ,RAC,$X, SY, SZ,EX ,EY ,EZ , 

ICX fCY tCZ v IBCOY ) 

C 

C THIS SUBROL TINE COMPETES THE PGSITICN GF THE SUN, EARTH, AND 

C CANCPUS IN THE PL AN ETOC 6NTR IC , PLANET EQUATCR, COORDINATE SYSTEM AND 

C WRITES CATA. THIS RCUTINE CALLS SUBROUTINES EEARTH, EMARS, EVENUS, PRECES 

C LATLNG, DOT, RECEC, RECVEC, AND REQMEQ. 

C 

C JD — JULIAN CATE AT TIME OF INTEREST 

C IEOLY - CONTROL INTEGER. 2 IMPLIES VENUS, A IMPLIES MARS. 

C CECS.RAS - CECLINAT ION ANC RIGHT ASCENSION CF THE SUN. 

C CECE,RAE - DECLINATION ANC RIGHT ASCENSION OF THE EARTH. 

C DEC C*RAC - CECLINATIGN AND RIGHT ASCENSION OF CANOPUS. 

C S X * S Y* S Z - UNIT VECTOR FRCM THE PLANET TO THE SUN. 

C EX , EY , EZ - LN IT VECTOR FRCM THE PLANET TO TEE EARTH. 

C CXtCYtCZ - IN IT VECTOR FRCM THE PLANET TO CANOPUS. 

C 

REAL JC 

RC=57 .295779513C823 
C 

CALL E EARTH! JD, >HE, YHE , ZHE , DXH E , DY h E ,DZ HE I 

CALL PRECES 12432 282 ., - .C EC34C592 , . 6C342 839 ,*.79513092 , JD ,C XE ,CYE ,C 
1ZE ) 

C 

IF € IB0DY.EQ.4) GO TC 2 
IF( IB0DY.E0.2 ) GO TC 1 

1 CALL EVENUS(JC,XHPtYHP,ZHP,DXHP,DYHP,DZHP) 

GO TO 3 

2 CALL EMARS1 JD, XHP,YFP,ZHP ,CXHP,DYHP,DZHP) 

C 

3 XHPE=XHE-XHP 
YHPE=YHE-YHF 
ZHP£=ZHE-ZHP 

RS£=$CRT<XFe**2*YFE<*2*ZHE**2) 

R SP = SORT ( XHP**2+YHF**2 + ZHP**2) 

R FE= SCRT 1XHFE** 2 + Y HP E** 2 *ZH PE **2 ) 

St X = XHE /R SE 
SEY=YHE/R SE 
SEZ=ZHE/RSE 
SP X= XHP /R SP 
SPY-Y HP/RSP 
SP Z- ZHP/R SP 

PE X- XHPE/RPE 
FEY=YHPE/RPE 
PEZ ^ZEPE/PPE 

CALL L ATL NG ( S EX , S E Y ,SEZ ,EHLAT , EHLONG ) 

CALL LATLNG (SFX , SPY , SPZ ,P HL AT, PHLONG I 
CALL COTC SEX,SEY, SEZ, SPX, SPY ,SPZ,ESPJ 
CALL DOT! SEX, SEY, SEZ , PE X, FEY, PEZ, SEP I 
CALL DCTC SPX,SPY,$P2,-PEX,-PEY,-PEZ,SPE) 

CALL RECEOl JD, — SPX, -SPY ,-SPZ ,SX£,SYE,SZE ) 

CALL RECEOl JC , P EX , P EY , PEZ , E XE, EY E , E ZE I 
C 
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I F ( JBODY.EQ .2 I GC 1C k 
IF t 1B0CY.EQ.9) GO TC 5 

k CALL RECVECl JD, I XE , ‘YE ,SZE ,SX ,SY ,SZ ,DECS,RAS) 

CALL REQVEC(JC,EXE, EYE.EZE.EX ,EY ,EZ,DECE, RAE) 

CALL RECVEC (JC,CXE ,CYE,CZE ,C X , C Y, CZ , DECC , RAC) 

GC TC 6 

5 CALL RECMEQ(JC»SXE«SYE*SZE*SX*SY»SZ»DECS»RAS) 

CALL RECMECUC, EXE , EY E , EZ E , EX , EY , EZ ,DECE , R A E ) 

CALL RLCMEQt JC.CXE.CYE ,CZE,CX ,CY ,CZ • DECC, RAC ) 

C 

6 XPS=SX*RSP 
YPS = SY*RSP 
ZPS=SZ*RSP 
X PE = E X*RP E 
YP E=E Y*RPE 
Z P E= E Z *RP E 

C 

IF ( IB0DY.E0.2 ) GC 1C 7 
I F ( IBODY .EQ • k ) GC TC 8 

7 WRITE ( 6 ,101 ) XHF ,YH F , ZHP .PHLAT , PHLONG, XHE , YFE , ZHE , EHL AT , EHLCNG 
WR IT El 6,102) XFS.YFi , Z P S , CE CS ,R AS , XPE , YPE , Z FE , DECE , R AE 

WRITE (6, 103) S X, S Y, S Z, EX, EY,EZ,CX, CY.CZ 
WRITE(6,1C4) FPE.RSE ,RSP ,5EP,ESP,SPE 
GL TC ECC 

8 WRITEI6 ,105 ) XF F ,YH P ,ZHF, FHLAT .PHLONG , XHE ,YFE , ZHE ,EHLAT, EHLCNG 
WRITE! 6,1 C6) X F S , YP S , ZPS , C ECS ,R AS , XPE, YPE, Z FE , DECE, RAE 

WR I T E ( 6 *1 0 7 ) SX,SY,SZ,EX,EY ,EZ,CX,CY,CZ 
WRITE (6 ,10 8) RPE.RSE ,R S P , SEP ,E SP , S PE 
C 

101 FOFMAT(*CFELICCENTPIC ECLIPTIC (MEAN EQUINOX OF DATE ) *// 2 k X , *X ( 
1KM)*,8X,*Y (KM)*,9X,*Z ( K N) * , 9 X, *L AT I T UOE* , 7X , *LON G ITUDE*/ / * VENU 
2S*,£X,5E16.7//* EARTH* ,5X, 5E16. 7) 

102 FORMAT (*GAPHPODIOCENTRIC (VENUS EQUATOR, VENUS EQUINOX )*//24X, *X 

1 (KM)*,8X,*Y (KM)*,9X,*Z (KM)*,6X,* DECLINATION RIGHT ASCENSION* 
2//* SUN* ,7X, 5E16.7//* EARTH* ,5X ,5 E16 .7) 

103 F C F M AT ( 14H C UNIT VECT0RS//12H SUN 3E16.7//12H EARTH 3E 

116.T//12H CANOPLS 3E16.7) 

109 FORMAT ( 29H0VENUS-EARTH DISTANCE (KM) = E16.8//29H SUN-EARTH DISTAN 
ICE (KM) = E16.8//29H SUN-VENUS OISTANCE (KM) = E16.8//29H VENU 

2 S- E AR T I — SUN ANGLE = E16.8//29H VENUS-SUN-EARTH ANGLE = E16 

3.8//29H E AR TH-VENUS-SUN ANGLE = E16.8) 

105 FORMAT(*CFEL ICCENTPIC ECLIPTIC (MEAN ECUINOX OF 0 AT E ) *// 29X , *X ( 
1KM ) * »8X , *Y ( KM ) * ,9X ,*Z (KM)*,9X,*LATIT UDE* , 7X ,*LONG ITUDE*/ / * MARS 
2*,6X,5E16.7//* EARTH* ,5X,SE16.7> 

106 FORM AT ( *0 AR EOC E MR IC (MAPS EQUATOR, MARS EQU INOX >*//29X,*X (KM)* 
1 , 8 X , * Y ( KM ) * ,9X ,* Z (KM)*,£X,* DECLINATIOM RIGHT ASCENSION*//* S 
2UN*,7X,5E16.7//* E ARTH* , 5X, 5 El 6.7 I 

107 FORM AT ( 14h0 UNIT VECTORS// 12H SUN 3E16.7//12H EARTH 3E 

116.7//12H CANOPUS 3E16.7) 

1C 8 FORM AT ( 28HCMAR S-EAR TH DISTANCE (KM) = E16.8//28H SUN-EARTH CISTANC 
IE (KM) = E16.8//28H SUN-MARS DISTANCE (KM) = E16.8//23H MARS-EA 
2RTH- SUN ANGLE = E16.8//28H MAR S- SUN-E AR TH ANGLE = E16.8//2 

38H E ARTH-MAPS-SUN ANGLE = E16.8) 

C 

ECC RETLRN 
END 
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INPUTS AND OUTPUTS 

The definitions of all input and output parameters are given in this appendix along 
with a sample input and the corresponding program output. 

Definition of Input 

Definitions of the program input follow: 


Program 

symbol 

Mathematical 

symbol 

Units 

Definition 

JD 

JD a 

day 

Julian date of arrival at planet 

SVDECE 

6 S,EE 

deg 

Declination of the incoming hyperbolic asymptote 
at the planet in the Earth equinox and equator 
coordinate system 

SVRAE 

X S,EE 

deg 

Right ascension of the incoming hyperbolic 
asymptote at the planet in the Earth equinox 
and equator coordinate system 

VINE 

Vco 

km/sec 

Hyperbolic excess velocity at the planet 

HA 

h a 

km 

Apoapsis altitude of the orbit about the planet 

HP 

h P 

km 

Periapsis altitude of the orbit about the planet 

BFRST 

ft 

deg 

The first orientation angle of the orbital plane 
to be considered where 0° % /3j ^ 360°. 

(See fig. 1.) If an iteration is performed on 
j3, the value of /3f must be an integer. A 
rational value may be input if only one incli- 
nation is considered. 

BLAST 

Vl 

deg 

The last orientation angle to be considered where 
0° s j3^ % 360° and is an integer; need not 

be input if only one inclination is considered. 

BSTEP 


deg 

Orientation step size where A/3 > 0; A/3 must 

not be input if only one inclination is considered. 

TLAST 

h 

day 

The last time in the planet orbit to be considered 
where t^ is measured from arrival; tj is 
considered zero if not input. 

TSTEP 

At 

day 

Time step size in the planet orbit; At need only 
be input if t^ is input. 

XJ20 

J 20 

none 

Second zonal harmonic of the planet; J 20 is 
considered zero if not input. (For Mars, 
J 20 = 0.002075.) 

IBODY 

IBODY 

none 

Integer indicator 

IBODY = 2 for Venus 
IBODY * 4 for Mars 

SUN1 

*1 

deg 

First Sun angle. The angle between the planet- Sun 
vector and the planet- spacecraft vector where 
0°<^; i// need only be input if Sun angle date 
is desired, in which case up to six different 
values of 1 // may be input. 

SUN 2 

4*2 

deg 

Second Sun angle 

SUN3 

^3 

deg 

Third Sun angle 

SUN4 


deg 

Fourth Sim angle 

SUNS 

^5 

deg 

Fifth Sun angle 

SUN 6 

^6 

deg 

Sixth Sun angle 


55 



APPENDIX C 


Sample Input 

The program input is loaded by using a FORTRAN IV namelist. A sample set of 
input data is as follows: 

$CASE JD = 2441533.5, 

SVDECE = 62.94, 

SVRAE = 120.12, 

VINF = 4.33, 

HA = 20000, 

HP = 1000, 

BFRST = 0, 

BLAST = 60, 

BSTEP = 10, 

IBODY = 2, 

SUN1 = 60, SUN2 = 70, SUN3 = 80, SUN4 = 900 
Definition of Output 

Definitions of the program output follow: 


Program 

symbol 

Mathematical 

symbol 

Units 

Definition 

Julian date 

JD C 

day 

Current Julian date 

Arrival date 



The calendar date of arrival at the planet 

SVDECE 

6 S,EE 

deg 

Defined under input (see SVDECE) 

SVRAE 

a s,ee 

deg 

Defined under input (see SVRAE) 

VINF 

Voo 

km/sec 

Defined under input (see VINF) 


The following output parameters are referenced to the mean planet equinox and 
equator of date coordinate system. They are listed in the same order in which they 
appear in the sample output. 


Program 

symbol 

Mathematical 

symbol 

Units 

Definition 

SVDECP 

5 s,pe 

deg 

Declination of the incoming hyperbolic 
asymptote at the planet 

SVRAP 

X S,PE 

deg 

Right ascension of the incoming hyper- 
bolic asymptote at the planet 
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Program 

symbol 

Mathematical 

symbol 

Units 

HA 

h a 

km 

HP 

h p 

km 

JD 

JD C 

day 

PLAT 

6 p 

deg 

PLONG 

X P 

deg 

PX 

p x 

none 

PY 

P Y 

none 

PZ 

P Z 

none 

QX 

«X 

none 

QY 

Qy 

none 

QZ 

Qz 

none 

WX 

w x 

none 

WY 

Wy 

none 

WZ 

W Z 

none 

VELgUN 

a 

deg 

SMA 

a 

km 

ECC 

e 

none 

INC 

i 

deg 

ARGPER 

O) 

deg 

ARGNOD 

ft 

deg 

DECSUN 

6 0 

deg 

RASUN 


deg 

DECETH 

6* 

deg 

RAETH 


deg 


Definition 

Defined under input (see HA) 

Defined under input (see HP) 

Current Julian date 

Latitude of periapsis 

Longitude of periapsis 

X- component of the P unit vector of 
the PQW triad 

Y- component 

Z- component 

X- component of the Q unit vector of 
the PQW triad 

Y - component 

Z- component 

X-component of the W unit vector of 
the PQW triad 

Y- component 

Z- component 

Angle between velocity vector at periap- 
sis and the Sun vector 

Semimajor axis of orbit about the planet 

Eccentricity 

Inclination 

Argument of periapsis 
Longitude of ascending node 
Declination of the Sun 
Right ascension of the Sun 
Declination of the Earth 
Right ascension of the Earth 
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Program 

symbol 

Mathematical 

symbol 

Units 

XSUN 

x © 

none 

YSUN 

y 0 

none 

ZSUN 

z o 

none 

XEARTH 

x® 

none 

YEARTH 

y® 

none 

ZEARTH 

z ® 

none 

PERIOD 

p 

hour 

XCANPS 

x c 

none 

YCANPS 

y c 

none 

ZCANPS 

z c 

none 

ZAP 

ZAP 

deg 

SVP 

<P 

deg 

VPE 

V pe 

km/sec 

VAE 

Vae 

km/sec 

VPH 

Vph 

km/sec 

DELV 

AV 

km/sec 

XJ20 

J 20 

none 

TIME 

t 

min 

T.A. 

f 

deg 


Definition 

X- component of the unit vector from the 
planet toward the Sun 

Y- component 

Z- component 

X-component of the unit vector from the 
planet toward the Earth 

Y- component 

Z- component 

Period of the orbit about the planet 

X-component of the unit vector from the 
planet toward the star Canopus 

Y- component 

Z- component 

Angle at the planet between the Sun vec- 
tor and the incoming hyperbolic 
asymptote (S-vector) 

Angle between incoming hyperbolic 
asymptote and the vector from the 
planet to periapsis 

Velocity at periapsis of the elliptical 
orbit about the planet 

Velocity at apoapsis of the elliptical 
orbit about the planet 

Velocity at periapsis on the hyperbola 

Deboost velocity (Vp^ - Vp e ) 

Defined under input (see XJ20) 

Time in orbit measured from periapsis 
to the point in question 

True anomaly of the point in question 
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Program 

symbol 

Mathematical 

symbol 

Units 

Definition 

ALT 

h 

km 

Altitude of the point in question 

DEC 

6 

deg 

Declination of the point in question 

RA 

X 

deg 

Right ascension of the point in question 

V/H 

V/h 

sec"l 

Ratio of horizontal velocity to altitude at 




point in question 
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Sample Output 

Samples of the output as received from the program follow: 


julian cm 

CALENCAR [HE 
ARRIVAL DATE 


2441523.5 

YEAR MONTH DAY 

72 6 4 

72 8 4 


HELIOCENTRIC ECLIPTIC (MEAN ECUINCX OF CATE) 

X (KM) Y (KM) l (KM) 

VENUS 1.0202678E+CE - 3. 4 5C 97 E2E + C7 -6 . 420 3809E +C6 

EARTH 1 .CC9E493E+C8 - 1. 12 298 69E +C8 0. 

APHRCCICCENTRIC (VENDS ECUATQR, VENUS ECU INOX ) 

X (KM) Y (KM) Z (KM) 

SUN -8.18C2037E+C7 7. 13 7269 1E + C7 7. 8121254E+C6 

EARTH -2.269E443E+C7 -7. 1997858E+07 -4. 5053729E+C5 

UNIT VECTORS 

SUN -7 • 5 1 56 <2 1E-C 1 

EARTH 
CANOFUS 

VENUS-EARTH DISTANCE (KM) 

SUN-EARTH CISTANCE (KM) 

SUN- V ENL S CISTANCE (KM) 

VENUS-EARTH-SUN' ANGLE 
VENLS-SLN-EfRTH ANCLE 
E ARTH-V EN l S-SUN ANGLE 


C.5574357E-C1 7 . 1 774666E-C2 

-4.135042 6E-01 -9. 1C ‘6 4 3 6E-G1 -5. 69749 1 1E-0 2 

8.1255C59E-C2 3.12C3492E-C 1 -9.4625934E-C 1 

7.9C764357E+C7 
1.51771 370E+C8 
1.08642379E+C8 
4.23967 36 3E+ 01 
2.994 40 75CE+01 
1.06659 189E + 02 


LATITUDE 

-3.3817195E+00 

0 . 


DECLINAT ION 
4*1 159245E+00 
-3. 2644396 E-01 


LONGITUDE 
1.8 518755E+01 
4 • 8288874 E+0 1 

RIGHT ASCENSION 
1 . 3889 514E+02 
1.1442558E+02 
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BET A 

50.CCCC 

(VENUS EQUATOR » VENUS EQUINOX) 

TIME FAST ARRIVAL = 

G.O 

DAYS 


SVCECE= 

6 .29400000 E+0 1 

SVRAE = 

1 .2C120000E+C& 

VINF = 

4.33000000E+00 

S V CECP = 

4.5751 5624E+01 

SVRAP = 

8.70988412E+01 

HA 

2.CCCC0C00E+04 

HP 

1. OOOOOOOOE + 03 

JO 

2 .441 53 350 E +06 

PLAT = 

1 .84920 24 8E + C1 

PLONG = 

4. 29 216827E+01 

PX 

6 • 83098535 E-0 1 

PY 

6 .5785 8492 E-0 1 

FZ 

3. 17172819 E- 01 

OX 

-6.2813579 2E -01 

QY 

3 • 264 3 1033 E -01 

QZ 

6.97298711E-01 

ux 

3.55188627E-01 

V*Y 

-6. 7872305 5E -01 

WZ 

6. 4278761 OE -01 

VELSUN= 

4. 19521438E+01 

SMA 

1 • 6585GCC CE *04 

ECC 

5. 7 280675 3E -01 

INC 

5. OOOOOOOOE +01 

ARGPER* 

2.44588570E+01 

APGNCO= 

2 .7623927 4E +0 1 

DECSUN= 

4.11 5 9245 2E + 00 

RASUN = 

1 • 38895138E +C2 

OECETH= 

-3 • 264439 62 E-0 1 

R4ETF = 

-1 .14425580E+Q2 

X SUN = 

- 7 • 51 5642 1 IE -01 

Y SUN = 

6. 55743573E-01 

ZS UN = 

7. 177466 56E-02 

XE APTH= 

-4.1350426CE-01 

YEA RTH = 

-9.1C46435 8E-0 1 

ZE ARTH= 

-5. 697491 14E-C 3 

PEP IOC= 

6. 540438 60E+00 

XCANFS = 

E.1255C588E-02 

YC ANP S = 

3 .13034915E-01 

ZC ANPS = 

-9 .46259344E -01 

ZAP 

6.1 1941 244E +01 

SVP 

4 .478392 2 7E + Cl 

VPE 

8 • 49 202 89 OE + 00 

VAE 

2. 3065372 7E+C0 

VPH 

1.0509 5498E+01 

DELV = 

2 • C 1 7 52C £ 5E + 0C 

> J20 = 

0 . 










TIME 

T . A . 

ALT 

DEC 

RA 

V/H 

SUN ANGLE 60.0 C EC AND 

DEC, 

S/C 

A SC 

12. 30 

46.64 

1912.76 

46.45 

89.58 

.00393299 

SUN ANGLE 60.0 DEG AND 

INCt 

S/C 

CSC 

91.32 

142. 47 

14334.23 

9.97 

-160.86 

.00020556 

SUN ANGLE 70 .C DEG AND 

DEC » 

S/C 

A SC 

8.03 

31.84 

1410.80 

39.59 

71.57 

.00568941 

SUN ANGLE 70.0 C EG AND 

INC, 

S/C 

CSC 

126.21 

157.27 

17539.44 

-1 .32 

-151.27 

.00014520 

SUN ANGLE 80.0 DEG AND 

DEC, 

S/C 

A SC 

4.42 

18.01 

1128.78 

31.15 

58.10 

.00738890 

SUN ANGLE 80.0 OEG AND 

INC, 

S/C 

CSC 

167.23 

171. 10 

19585.12 

-11 .85 

-142.23 

.00011967 

SUN ANGLE 90.0 DEG AND 

CEC, 

S/C 

ASC 

1.11 

4.55 

1C08.16 

21.81 

47.24 

.00841363 

SUN ANGLE 90.0 DEG AND 

INC, 

S/C 

CSC 

-181. 27 

-175.45 

19890.03 

-21.81 

-132.76 

.00011646 

SUN CCCULT A TION TIME 




0.00 






ENTER SUN GCCULTA7 ICN 




0.00 

C. CC 

0.00 

0.00 

0.00 


EXIT SUN CCCULTATICN 




0.00 

0.00 

0.00 

0.00 

0.00 


EARTH GCCUL T AT ION TIME 




22.97 






ENTER EARTH CCCULTATION 




-10.97 

317.75 

1740.23 

-13.53 

15.97 


EXIT EARTH OCCULTAT ION 




11.99 

45.63 

1871.47 

46.08 

88.23 


CANCFLS CCCULTATICN TIME 



0.00 






ENTER CANCPUS OCCULTATION 



0.00 

0.00 

0.00 

0.00 

0.00 


EXIT CANCPLS CCCULTATICN 



0.00 

C.00 

0.00 

0.00 

0.00 
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CT> 

CO 



BETA 

= 

60.0000 (VENUS 

EGUATCR, 

VENUS EQUINOX) 




TIME PAST ARRIVAL = 

C.O DAYS 





SVC ECE = 

6 .294C0CCCE+01 

SVRAE = 

1.2C12000CE+02 

VINF = 

4. 33000000E + 0C 

S VCECP= 

4 #5751 56 24E +01 

SVRAP = 

£ • 7CS88 4 12E + 01 

HA 

2.CCCC0000E+04 

HP 

1 .OOOOOOOOE +03 

JC 

2 .441 5 3350 E +06 

PLAT * 

9 .53023C77E +00 

PLDNG = 

5.621239 76E +01 

PX 

5* 47009 130E -Cl 

PY 

8. 20590189 E-01 

91 

] • 6 55679 7 3E-C 1 

GX 

-5.01C16087E-01 

QY 

I • 62467 333E-0 1 

QZ 

8.50051320E-01 

l*X 

6.7C644386E-01 

WY 

-5.47938C51E -01 

WZ 

5 • CCOOOOOOE-O 1 

VE LSUN = 

5*703717 16E+01 

SMA 

1 .65E5CCC0E+O4 

ECC 

5.728C6753E-01 

INC 

6. OOOOOOOOE +01 

ARGPER- 

1 • 10 217390 E+0 1 

ARGNC C= 

5.C75Q0715E+01 

□ EC SUN = 

4. 1 159 2452E+00 

RASUN = 

1 • 3 8895 138 E +0 2 

DECET H= 

-3.26443962E-0 1 

RAETH ^ 

-1.14425580E+02 

XSUN = 

-7. 515 64 211 E- 01 

Y SUN = 

6. 55 743 573 E -0 1 

ZS UN = 

7 • 1 77466 56E-02 

XE ARTH= 

“*•135 04 26 CE -01 

Y EA RTH = 

— 9 . 1 C4 E4 35 8E-01 

zearth= 

-5.697 49 1 1 4E —0 3 

PERIOD= 

6.540438 60E+00 

XC AN FS = 

f . 125505E6E-C2 

YC ANP S = 

3* 1303491 5E -01 

ZC ANP 5- 

-9*4625 9344 E -0 I 

ZAP 

6*119412 44E +01 

SVP 

4 .4783932 7E + 0 1 

VP E 

8.492 02890E+00 

VAE 

2. 3065372 7E + 00 

VP H 

1.05095498E+01 

CELV = 

2 .01752C85E+00 

X J 20 = 

0 . 






SUN 

ANGLE 

6 C 

.0 

DEC 

AND 

INC , 

S/C 

CSC 

SUN 

ANGL E 

60 

.0 

CEG 

AND 

DEC, 

S/C 

A sc 

SUN 

ANGLE 

70 

.0 

DEG 

AND 

INC, 

S/C 

CSC 

SUN 

ANGL E 

70 

• 0 

DEG 

AND 

DEC , 

S/C 

ASC 

SUN 

ANGLE 

80 

.0 

CEG 

AND 

INC, 

S/C 

DSC 

SUN 

ANGLE 

80 

*0 

DEG 

AND 

CEC, 

S/C 

ASC 

SUN 

ANGLE 

90 

.0 

DEG 

AND 

INC, 

S/C 

DSC 

SUN 

ancl e 

90 

.0 

DEG 

AND 

CEC, 

S/C 

ASC 


SUN CCCUIT #T ICN TINE 
ENTER SUN CCCULTAT ION 
EXIT SUN CCCULTATICN 


EARTE OC CULT AT ION TINE 
ENTER EARTH OCCULT AT ION 
EXIT E ART N OCCULTATION 

CANOPUS CCCULTATICN TIME 
ENTER CANOPUS OCCULTATION 
EXIT CANCELS OCCULTATION 


TIME 

T . A. 

ALT 

DEC 

RA 

V/H 

35.46 

102. 76 

6671.89 

52.42 

-177.86 

.00070690 

12.92 

48.61 

1997.22 

48 .35 

91.22 

.00 37 27 31 

66.78 

128. 16 

11162.31 

34.48 

-152.61 

.00031252 

5.75 

23.21 

1215.08 

29.15 

69.54 

• 00678292 

102.41 

147.67 

15510. 85 

18 .34 

-140.29 

.00017962 

.90 

3.70 

1005.37 

12.71 

58.23 

.00844028 

1 50. 38 

165. 68 

18957.04 

2.85 

-130.90 

.00012674 

-3.50 

-14.32 

1081.06 

-2 .85 

49.10 

.00776639 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 


0 .00 

C.00 

0.00 

0.00 

0.00 


25.38 

-13. 42 

3 0 9 . 81 

2063.20 

-33. 16 

28.59 


11.96 

45. 52 

1867.02 

46.2 6 

87. 86 


0.00 

0.00 

0.00 

0.00 

0 .00 

0.00 


0.00 

0.00 

0.00 

0.00 

0.00 
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